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SYNOPSIS 


This study analyzes the interaction between shock wave 
and laminar boundary layer on a flat plate numerically, mainly 
emphasizing the influence of suction on the separated region. 

The complete compressible Navier-Stokes equations are 
applied to the problem in two parallel finite difference formu- 
lations, by Brailovskaya and MacCormack. 

In the present scheme the outer edge of the computational 
region is placed just outside the viscous layer; assuming a 
simple wave type of the flow there, the outer boundary condition 
is computed by extension of the flow field characteristics. 

This allows the reduction of the computational region, saving 
computer time and storage. 

To shorten computation time even more, an approximate 
shape of the boundary layer form along the wall, having a 
small separation bubble below the shock, is included in the 
Initial conditions. 

The computational field is 150 & Q long and 3 <$ 0 wide (S Q - 

boundary layer thickness at the entry) consisting of 76x25 mesh 
nodes. The computer time required for obtaining a converged 
accurate solution (reached within about 200 time steps) Is of 
the order or 4 CPU minutes on the IBM 370/168 computer. Good 
agreement with known results is achieved. 

This relatively rapid solution, differing only slightly 
from the original MacCormack results that converge more slowly, 
as in other previous studies, allows the practical application 
of this scheme to many cases of different flow conditions 
within reasonable computation time. 

In developing the computational scheme, several auxili- 
ary methods were tested (though not Included in the final 
scheme), as follows: artificial viscosity, non-uniform grid 
meshes, fourth order differences, exclusion of the shock wave 
from the computational field, and variation of initial and 
boundary conditions. 

The present Brailovskaya scheme has been used for param- 
etrical study of the interaction flT3w field for the following 
range of parameters: entry Mach number 2 < M Q < 4.5; downstream 
to upstream pressure ratio 1.2 5. P e /P 0 — 3-2; and impingement 
Reynolds number 104 < Re x _ < 10°. The influence of thermal 

conductivity (cooling and warming) and mass transfer (suction 
and injection) along the wall has been tested as well. 
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The results analyzed primarily were the size of the 
separated region, wall pressure distribution and boundary layer 
profiles ^ 

The applicable interpretation of the separated region is 
very important since it greatly affects the performance 
efficiency of aircraft components, such as engine inlets and 
control surfaces, where such an interaction takes place. 

Wall suction has been found to be a useful means for 
reducing separation until it disappears completely. A compre- 
hensive analysis was made for evaluation of the amount of 
suction 'needed- to prevent boundary layer separation. Given 
the flow conditions previously mentioned, the required suction 
lies in the range between 0 > > -0.04 (i.e. up to 4 % of 

the free stream velocity). 
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NUMERICAL SOLUTION FOR THE INTERACTION OF SHOCK WAVE 
WITH LAMINAR BOUNDARY LAYER IN TWO-DIMENSIONAL FLOW 

ON A FLAT PLATE 


Uriel Landau 




1. Introduction /10* 

1.1. General 

This research attempts to analyze the phenomenon of inter- 
action between a shock wave entering the laminar boundary layer 
of a flat plate, as well as some of the accompanying practical 
projections, through the employment of the advanced methods of 
computation available today in the field of numerical calcula- 
tion of flows. 

From a practical point of view this interaction phenomenon 
(in which shock waves, expansion and compression fans, the 
boundary layer and the flow separation zone are involved) occurs 
at the inputs of jet engines and at the control surfaces of 
supersonic aircraft (when there is tubbulent flow, as is true 
in some cases). This interaction has considerable influence 
on the flight performance and therefore must be considered in 
the design of components that are directly or indirectly in- 
volved with it so as to improve that performance during all 
flight conditions as much as possible. Investigations of the , 

Influence of interaction have, until recently, been carried out ^ 
only experimentally (in a flight tunnel or in actual flight) 
because of the absence of any reliable and accurate analytical 
methods, which were, moreover, considered too expensive and 
too drawn out. 

Recently, with the rapid development of computers, it has 
become possible gradually to provide numerical solutions to 
more and more problems involving voluminous computation. The 
time has thus arrived to develop a method of efficient and 
relatively inexpensive solutions for problems of this nature. 

The numerical solutions that have appeared in recent years 
on this subject constitute an important step forward and 
while much remains to be done to Improve and expand their 
practical application, they represent a good basis from which 
to continue the task. 


This study recommends a method of computation that is 
based on the expansion and utilization of several presently 
known methods and which is basically a numerical solution of 
the complete (Navier-Stokes equations) equations for the 
field of flow close to the wall, including all its zones, 
as a single unit. By accepting an accuracy that falls only 
a little short of what is presently available, a convergence 
can be achieved after a relatively short period of computation 

^Numbers in the margin indicate pagination in the foreign text. 
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(a few minutes on the CPU) and a practical method is realized 
through which the computation method can be applied to a 
large scope of different flow conditions (which have never 
before been examined through any similar type of computation) 
with two primary aims: 

(A) Examination of the influence of various parameters 
on the nature of interaction; 

(B) Investigation of practical ways (especially suction) 
to prevent separation of the boundary layer during the inter- 
action. The results of the calculations show that there is 
room for the continuation of the development of methods of 
this kind and for their expansion to additional cases (such as 
flows around corners and steps and the entry of a turbulence 
model into the calculations). 

1.2. Qualitative Description of the Phenomenon 

The phenomenon of interaction, which is shown schematically 
in Fig. 1, occurs as the outcome of a shock wave impinging at 
a slant angle on the laminar boundary layer* of a flat plate.** 
The shock wave intersects the boundary layer while going 
through a process of thickening and bending and reaches the 
line where the flow is sonic. The pressure increase that 
follows the shock wave is transferred backward, upstream (up- 
stream influence) through the subsonic portion of the boundary 
layer. This increase in pressure causes a gradual increase in 
the thickness of the layer just before the point of entrance 
of the shock wave. Where the shock wave is strong enough this 
may lead to separation of the boundary layer at the point where 
the internal friction at the plate becomes zero. At the point 
of separation the separation zone starts, which is detached 
from the boundary layer by the flow line and in which a 
recirculating flow is maintained. 

Circulation of the flow in the area before the entrance 
of the shock wave creates a spread of the compression waves 
(because of the change in boundary layer thickness) which 
unite into a shock wave of separation at considerable distance 
from the plate. 

The separation zone causes the return of the wave spread- 
ing since this area cannot sustain a sudden pressure increase. 
Most of the pressure rise generated by the impinging shock wave 


*A similar phenomenon is accepted in the turbulent boundary 
layer but the changes in this field of flow are smaller. 

**The above-mentioned interaction occurs at sufficient distance 
from the leading edge of the plate to permit the disregard of 
the influence of interaction between the boundary layer and the 
shock wave that emanates from the leading edge. 
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is attenuated by the fanning out of the expansion waves. The 
incident . wave and the spreading of the expansion waves return 
the direction of the external flow toward the plate. The 
turn of the flow from this condition until it again becomes 
parallel to the plate, at the downward slope of the flow, 
causes the fanning out of the compression waves, which unite 
far from the flat plate into the reattachment shock wave. 

Two shock waves generated by the compression fans unite further 
still into one wave, the reflected shock wave. 

The turn of the flow toward the plate reduces the thickness 712 
of the separation zone gradually up to the point where the 
boundary layer is again attached to the plate, while its thick- 
nes continues to decrease to its minimum in the zone where the 
external flow is parallel to the flat plate. During the down- 
ward slope of the flow from this zone the boundary layer increases 
again. 

All the above-mentioned phenomena depend primarily on the 
strength of the incident shock wave with regard to the Mach and 
Reynolds numbers of the flow in the interaction zone and to the 
laminar or turbulent nature of the boundary layer. 

The order of magnitude of the interaction zone in which 
the majority of the above-mentioned changes take place is a 
7"" length of 100 S 0 and a width of 3 S 0 (6 0 being the width of the 

boundary layer at the input cross section). Outside this area 
the influence of the interaction on the external field of flow 
and on the boundary layer is negligible. 

From the point of view of the mathematical nature of the 
flow equations, there are three regions in the interaction zone: 

— the parabolic region: the boundary layer (without 
separation zone and shock wave vicinity); 

— the hyperbolic region: the external flow outside the 
boundary layer; 

— the elliptical region: the separation zone and the 
vicinity of the shock wave (along its incidence to the 
boundary layer). 

The last two regions do not conform to the boundary layer, 
mainly because they contain strong gradients of the same 
order of magnitude both in the direction of the flow and at 
right angles to it. 

To avoid the use of three sets of separate equations for 
each region and the stitching together of their solutions for 
N the lines of contact between them, the complete Navier-Stokes 
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equations are employed which are applicable to all regions. 

This set .of equations in its non-continuous form is parabolic 
and becomes elliptical in its continuous form (the asymptotic 
solution is found from the continuous form). This arrangement 
requires initial conditions for each field and boundary condi- 
tions for all boundaries of the field of flow (since the problem 
being solved is an Initial Value - Boundary Value problem). 

Each field of interaction constitutes thus a single unit with 
all its regions, at least for the purpose of computation. 


I.3» General Background of the Research* 
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The subject of interaction between the shock wave and the 
laminar boundary layer at a flat plate for two-dimensional flow 
was first investigated experimentally (starting at the end of 
the 19^0' s) and afterwards through various analytical methods 
which gradually improved. 

The first analytical solutions, which used the method of 
integrating the momentum along the width of the boundary layer, 
produced results (applicable only to special cases where experi- 
mental results already existed) that depended largely on the 
choice of empirical coefficients and on various mathematical 
assumptions . 

Much progress was made by means of a solution of the boundary 
layer equation through finite differences, which represented a 
more general solution, included fewer assumptions, and provided 
results that were better but that still depended on the close 
formulation of the interaction mechanism; there also still 
remained the requirement to adapt the solutions to the various 
zones of flow. 

A new direction for the improvement of the computation was 
taken in recent years with the expanding use of computers, which 
led to the founding of the new field of Computational Fluid 
Dynamics. Within the realm of this field many flow problems 
are being solved by means of numerical solution of sets of 
equations and boundary conditions that represent the physical 
phenomena far more accurately. Thus the approach developed, 
which is based on the solution of the entire field of flow 
(i.e. boundary layer, external flow, the vicinity of the inci- 
dent shock wave and the separation zone) through finite 
differences integration of the complete equations (continuity 
equations, the momentum in the x and y directions and the 
energy) with suitable boundary conditions. This system is 
solved through an iterative process (by integration of the 
equations vs. time to obtain an asymptotic solution). 

A few of those studies that used this method achieved 
results that were better than previous ones, the computation 


*The subject is treated at length and with references in Part 2. 
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became more general and easier to perform. Even though there 
is room for improvement and perfection of the computation 
method, there is no doubt that this approach provides far more 
trustworthy solutions for physical phenomena like the above- 
mentioned interaction. 

The main subjects still open for additional improvement 
(some of which was partly accomplished as one of the goals of 
this research) are: 

a) Perfection of the method of computation to increase the 
accuracy of results and reduce the field of computation to the /l4 
immediate vicinity of the interaction zone, as well as to 

speed up the rate of convergences (reduction of computing time). 

b) Expansion of the range of computation to stronger inter- 
actions than had been computed in the past and setting up of 
parametric research oh the influence of various factors on the 
interaction. 

c) More accurate representation of the shock wave (in the 
present methods its width is several times that of its physical 
width, which leads to some distortion of results). 

d) Examination of actual projections of the phenomenon, 

the chief one of which is to find a way for reducing the separa- 
tion zone (especially through suction), which is responsible for 
the inefficient performance of the parts of an aircraft on 
which this interaction takes place. 

These subjects, whose investigation represents the aims of 
this and other research studies presently carried on in this 
field, have been dealt with in three principal stages: 

a) A definition of the differences integration, of the 
computation method for the equations and of the boundary condi- 
tions, that provides a resounding solution to the problem; 
various approaches and techniques, from which optimum combina- 
tions are formed for the computation method. Among the 
approaches that have been examined were: "continuous" and "non- 
continuous" treatment of the shock wave (in the "non-continuous" 
treatment the shock wave is removed from the field of computa- 
tion) , employment of artificial viscosity, computation of 
differences of the second order and of the fourth order, a 
uniform computation grid and a non-uniform one, varying 
boundary conditions and initial conditions. By means of these 
approaches two final schemes of computation were prepared, 
which are based on the methods of Brailovskaya and MacCormack. 

b) Application of the method of computation to a series 
of variable flow conditions, for the investigation of the 
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characteristics of interaction and the factors that influence 
it, when .the principal results in question are the magnitude 
of the separated flow zone, the shape of division of pressures 
and the friction along the length of the plate, as well as 
the boundary layer profiles. 

c) Detailed examination of the influences of suction on 
the separation zone under various flow conditions, with the aim 
of finding the required condition for diminution or total 
elimination of the separation zone. 

A detailed summary of the changes and innovations made 
during this study, with reference to existing studies, is pre- 
sented in Section 2.3 (after a survey of the references). 


n 


2. Background and Reference Survey /15 

The subject of interaction between shock waves and boundary 
layers and the methods for its solution is a voluminous one 
and branches off into many directions. Prom the wealth of 
references on the various aspects of the subject only a few 
of the principal ones were chosen, which are connected with 
the computation methods used and the results achieved in this 
research study. 

Later on two central subjects will be dealt with separately 
(each of them is subdivided into many secondary subjects): 

a) The subject of interaction between the shock wave and 
the laminar boundary layer of the flat plate, as investigated 
up to the present both experimentally and analytically. 

b) Various numerical methods for the solution of Navier- 
Stokes equations in compressive and viscous flow, which serve 
as the principal tool for numerical solution of problems of 
fluid dynamics. 

2.1. Interaction Between Shock Wave and Laminar Boundary Layer 
2.1.1. Background References 

The references described below are mostly survey papers 
that deal with the subject from various points of view and 
in its various stages of research development. 

The phenomenon of interaction, against the broad background 
of various flow phenomena, has been treated in a number of 
basic texts; Shapiro [93] deals with the subject in a broad 
analytical discussion of compressive flow with shock waves. 
Schllchting [92] (1968), [text illegible] and Stewartson [95] 
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(1964) treat the interaction phenomenon in a discussion of the 
characteristics of the boundary layer during compressive flow. 

Dorrance [23] (1972) deals with viscous hypersonic flow 
and ascribes the various projections of the phenomenon to 
chemical reactions of the fluid and/or the surface area. In 
all of these references a general and qualitative description 
of the phenomenon is given. 

A broad survey, accompanied by quantitative results from 
practical problems that Involve interactions between shock waves 
and the boundary layer on aircraft, is presented in the collec- 
tion of papers edited by Lachman [55] (1961) and in the study 
by Korkegi [53] (197D, in which the decreased efficiency of 
control surfaces and engine inlets on jet engines is described 
as due to the flow separation caused by interaction and which 
suggests ways for improvement of their performance. 

2.1.2. Experimental Studies /l6 

The early work dealing with interaction was experimental. 

This was so primarily because the physical mechanism of the 
phenomenon was not well enough understood at first and the 
analytical and computational tools were lacking for accurate 
treatment of the subject. 

One of the earliest studies was conducted by Liepman et al. 
[6l] (1949) and offers a qualitative analysis of the phenomenon 
(based on Schlieren photographs) as well as a model of the flow 
in the interaction zone, by means of which the magnitude of 
the zone that is subject to back pressure influence was set 
(about 50 <$ Q ). Two additional investigations were by Barry 
et al. [4] (1950) and Gadd et al. [28] (1954), which also 
describe qualitative results and analyze the interaction mechanism 
for the first time. 

The following investigations already describe quantitative 
results for the pressure divisions, internal friction and thermal 
conductivity along the plate, for various Mach and Reynolds 
numbers and shock wave forces. Care is taken to maintain the 
two-dimensional and laminar characteristics of the flow under 
experimental conditions. 

Chapman et al. [l4] (1958) and Hakkinen et al [4l] (1959) 
provide results for average Mach numbers and shock wave forces 
(1.2 < M 0 < 2.5; 1.2 < P e /Po < 2.4) for interactions that 
occur when a shock wave is incident on the boundary layer of 
a flat plate. Needham, Stolery [75] (1966) and Lewis et al. 

[60] ( 1968 ) expanded the range of flow conditions (to M 0 = 10 
and p e /p Q =4) when the interaction occurs around a compression 
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corner and part of the flow (only for the strongest interactions) 
is turbulent. Those results are also of significance for this 
research since such an interaction is similar in its appearance 
and in its qualitative (as well as to some extent in its quanti- 
tative) results to the one caused by a shock wave impinging on 
the boundary layer of a flat plate. 

2.1.3. Analytical Studies 
2. 1.3.1. General Approach 

The first attempts at solution of the interaction were based 
on various approximations of the boundary layer in the vicinity 
of the zone where a weak shock wave is incident. In those 
first investigations by Lighthill [62] (1950), Lighthill [ 63 ] 

(1953) and Stewartson [95] (1955), the first analysis of the 
interaction mechanism and of the back pressure influence is made, 
through experiments with the building of an 'analytical model to /17 
explain the phenomenon and by using omissions and approximations. 

The results obtained were not accurate and depended too much on 
the approximations, while the method of computation could not 
be extended to more general cases. 

At a later date the approach (which remained up to date 
until the end of the 1960’s) to the solution of interactions 
through boundary layer equations, with the addition of an 
expression that tied the external flow to the boundary layer, 
was developed. Those solutions were mostly obtained through 
the integral moment method, though some employed the method of 
finite differences of the boundary layer equations. 

Solutions obtained through the integral moment method 
improved with time thanks to the perfection of computation tech- 
niques (such as different solutions for different zones and 
their integration, the use of semi-empirical coefficients, etc.), 
but their accuracy remained limited to low Mach numbers and weak 
shock waves; the solutions were also too dependent on the assump- 
tions of the method of computation and on the chosen coefficients. 

The method of finite differences for the solution of boundary 
layer equations perfected the computation and made it more general. 
Results improved and the range of solutions was broadened. But 
here, too, the disadvantage of using basic assumptions about 
the boundary layer In zones where they did not apply (the sepa- 
ration zone and the vicinity of the shock wave inside the 
boundary layer) was conspicuous. 

Lately, with the great development in computers and the 
investigation of numerical methods, a new approach has taken 
shape (which is also typical for other solutions in compressible 
flow) in which complete Navier-Stokes equations are solved by 
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finite differences, for the entire field of flow and its 
various zones. The solution becomes a more general one and is 
only influenced by the choice of flow parameters and by the 
boundary conditions of the problem. This approach, which was 
started under the auspices of a new field called Computational 
Fluid Dynamics, has already yielded excellent solutions for 
the interaction in a number of special cases. 

2. 1.3. 2. Solution of the Field of Flow Using Boundary 
Layer Equations 

The different methods available for this approach will be 
described briefly further on. A thorough survey of the main 
methods applicable for the solution of interaction through use 
of the boundary layer equations provides the following refer- 
ences: Holden [45] (1965); Morduchov [73] (1955); Holder, 

Gadd [46] (1955); Murphy [74] ( 1969 ); Georgeff [32] (1972); 

Charwat [15] (1970); Panov, Schvets [79] (1967). The references 
apply mainly to the integral moment method, and only Murphy [74] 
contains finite difference methods as well. The four last- 
named papers describe (through analysis of the results and 
their comparison with experimental results) the inaccuracy of 
the approach, which is due mainly to omission of the gradients 
in the y direction (according to the boundary layer equations) /l8 
in flow zones where it is not applicable (such as in the separa- 
tion zone and the zone where the shock wave is incident on 
the boundary layer). 

What the studies that share this approach have in common 
is their reliance on the concept of free interaction by Chap- 
man. According to it, the increase in pressure via the separa- 
tion zone up to the plateau of steady pressure does not depend 
on a mechanism that generates the interruption (like an 
entering shock wave, a forward step, a compression corner, etc.). 

Now we will take up the methods mentioned above, each one 
separately : 

2. 1.3. 2.1. The Integral Moment Method 

This approach is based on an assumption of the existence 
of boundary layer equations in which functions or tables for 
the dependence of unknowns in the y direction are posited, 
by means of a dependent parameter x; all of this based on 
experimental results or theoretical assumptions. (It becomes 
clear then that the accuracy of those functions or tables 
determines the accuracy of the computation result’s.) The 
resulting equations are multiplied by weighting functions of 
one of the variables (velocity for the Moment Equation and 
temperature for the Energy Equation) and are integrated in 
the y direction (along the width of the boundary layer) . The 
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result of this integration is a set of differential equations 
of the first order for the variables in the boundary layer 
that are functions of x. To this set is added an equation of 
the form 



which expresses the interaction between the pressure change due 
to the impingement of the shock wave and the change in the 
shape of the boundary layer. 

The resulting equations are integrated in the x direction 
as far as the decline of the flow where the adaptive equations 
become valid 
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The investigations by Glick [35] (I960) and Bray et al. 

[10] (1961) are based on the Integral Moment Scheme of Crocco, 
Lees [19] (1952), which deals with the interaction between 
dissipative flow and nearly isentropic flow. In those studies 
the continuous equations for the momentum and of the moment 
of the momentum, for the flow at an adiabatic wall, were con- 
sidered and separate velocity profiles were used in them, 
both for an attached and a separated boundary layer. The start 
of interaction was generated by disturbance through pressure. 

The influences of various parameters (like suction and cooling 
of the wall) on the solution, which is expressed by the pres- 
sure division, the friction and the form of the boundary 
layer’s development, were discussed. 

This line of thought was extended to the non-adiabatic 
flow (by inserting the enthalpy profile into the energy equation) 
by Klineberg, Lees [52] ( 1969 ); Holden [45] ( 1965 ); Lees, 

Reeves [59] (1964); Georgeff [32] (1974) and also by Horton 
[47] (1971) who discussed the aximetric problem. 

The investigations by Nielsen et al. [ 76 ] (1965); Goodwin 
et al. [ 38 ] (1967); Makofski [70] ( 1963 ) also deal with non- 
adiabatic flow, but they offer separate solutions for the 
separated zone and the boundary layer on both sides of the 
line of flow that separates them. In the third study the 



velocity profile is defined as a function of the friction and 
of the pressure gradient along the plate. In Blankenship's 
[6] (1967) paper the solution is obtained with the aid of 
small disturbances to generate interaction in tubes. 

The solutions demonstrated in the above studies are partly 
dedicated to problems of the external shock wave, which 
impinges on the boundary layer, and partly to problems of the 
pressure corner, two problems that are similar in their nature 
and in the results from their solutions. 

2. 1.3. 2. 2. The Finite Differences Method (Boundary 
Layer Equations^ 

This approach was developed later, thanks to developments 
in the use of computers. In the computation the boundary 
layer equations, together with the equation of interaction 
(2-1), constitute a series of parabolic equations that are 
solved numerically by an implicit method of finite differences. 
Similarly, the appropriate boundary conditions are formulated 
for upstream flow (the profile of initial flow) and for the 
length of the plate and the external flow as well. A solution 
through this method was offered by Reyhner, Flilge-Lotz [84] 
(1968). For the sake of including the separated (aliphatic) 
region in the range of computation of the (parabolic) equations, 
several omissions were made that do not influence the result 
much but which prevent unsteadiness in the process of numeri- 
cal computation. Except for that, no further empirical 
assumptions are made and the results are much better than those 
from the Integral Moment Method (when compared to experimental 
results) . 

A solution through finite differences, by Baum [ 5 ] (1968), 
must be mentioned, which relates to interaction by means of 
a rounded corner and in which the boundary layer equations 
contain a parameter for the surface curvature. 

2. 1.3-2. 3. Other Methods of Solution 

In addition to the methods referred to in the two previous 
paragraphs, various partial solutions were obtained (by means 
of approximations of equations), using other methods than 
those already mentioned: 


2. 1.3-2. 3.1. Division of the Boundary Layer into 720 

Subregions 

Based on Lighthill's [ 63 ] (1950) approach, the method of 
solution through division of the boundary layer into two 
sublayers and an approximation for the solution between them 
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was developed. Goodman [36] (195*0 suggested such a solution 
path In which Interaction is expressed by the insertion of a 
semi-empirical function for the pressure division around the 
region where the shock wave impinges. Ray [83] (1962) devel- 
oped semi-empirical formulas for the description of various 
results like pressure division, boundary layer thickness, 
etc., in this way. 

Rose [88] (1969) and Miller [72] (1973) perfected the 
method and solved the inner (viscous) part of the’ 'boundary 
layer through boundary layer equations and the outer (inviscid 
and rotational) part through a method of characteristics, with 
the two solutions matched to each other at the line of contact 
(which is the sonic line on which the shock wave impinges). 

For Rose's [88] solution an empirical value is needed for the 
length of the region that is influenced by back pressure. 

Brilliant, Adamson [11] (197*0 solved interaction in 
transonic flow through division of the field into several 
regions: the inner boundary layer through boundary layer equa- 
tions, the outer boundary layer and also the external flow 
through equations of the small separate disturbances. The 
three regions were matched up at their lines of contact. 

2. 1.3. 2. 3. 2. Various Semi-Empirical Methods 

In those methods various semi-empirical functions are 
employed and coefficients are computed by approximation methods 
to get specific results that characterize the interaction. 

That is how Erdos, Pallone [25] (1962); Popinski [82] 

(1965) and Gai [29] (1970), as well as Byrkin [12], did their 
work. Byrkin [12] obtained analytical approximate solutions 
for a number of special cases. Gai [29], who relied on compu- 
tation of the sublayer of the subsonic boundary, was mostly 
concerned with computing the distance over which the influence 
of interaction was felt upstream of it, for various Reynolds 
numbers . 

Still to be noted is Ackroyd's [1] (1969) work, which 
discussed analysis of the influence of various parameters on 
the mechanism of interaction. 

2. 1.3. 3- Solution of the Whole Flow Field Using /21 

Complete Navier Stokes Equations 

Thanks to the development of fast computers in recent years, 
it has become possible to solve the complete Navier-Stokes 
equations numerically. Many investigations of a wide range of 
subjects have recently been provided with solutions in this way. 
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solutions of the interaction between a shock wave impinging on 
the boundary layer have been obtained by Skoglund, Gay [9*0 
(1969) and MacCormack [ 67 ] (1971). Since those are the most 
advanced studies in this field and very close to the subject 
of this research, their discussion will be more detailed 
than that of preceding work. 

The two studies are based on developments in the explicit, 
two-stage difference schemes by Lax-Wendroff*, which helped 
improve convergence and increase accuracy especially near the 
wall. Skoglund, Gay [9*0 did a logarithmic transformation 
of the coordinates, which reduces the size of grid meshes in 
the vicinity of the incident shock wave on the boundary layer 
and increases their size gradually toward external flow in 
the y direction and toward upstream flow in the x direction. 
MacCormack [ 67 ] divided the grid into two regions; in the 
region near the plate the grid was divided into small spaces 
in the y direction and the distant region was divided into 
larger spaces (the meshes being uniform in each respective grid). 
A saving in computer time for the region close to the plate was 
achieved through subdivision of the computation by difference 
equations into x and y directions separately by a special 
technique. Similarly, the differences in changing directions 
(alternately forward and backward) were .computed. 

Skoglund, Gay [9**] added terms for artificial viscosity 
for the sake of stabilizing the solution in separated flow, 
while MacCormack [ 67 ] did not require it since in his differ- 
ence schemes there is already enough (numerical) artificial 
viscosity for stabilization of the solution. 

In both investigations the boundary conditions were 
formulated in a similar manner. Skoglund and Gay [9*0 
assumed an initial boundary layer profile per Chapman-Rubenstein 
at the entrance cross section, a linear extrapolation of the 
variables at the exit cross section, and along the plate no 
components of velocity, adiabatic wall temperature, as well 
as density from one of the equations of motion. MacCormack 
[ 67 ] assumed uniform flow conditions along the entire width 
of the entrance cross section (which, for him, was in front 
of the leading edge of the plate), zero gradients at the exit 
cross section and along the length of the plate no velocity 
components and zero gradients for pressure, temperature and 
density in the y direction. 

In both studies the external flow boundary was assumed to 
be far enough from the plate so that the reflected flow lines 
of the incident shock wave will exit via the exit cross section 


*See Paragraph 2.2.2. 
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and not via the plate. Conditions at the outer boundary were 
therefore determined as fixed according to the formulas for 
shock wave transition (from both sides of the incident shock 
wave). The shock wave itself was determined as having a 
width of 1-2 grid meshes in its initial form. 

In both investigations the numerical results were calculated 
for M Q = 2; Re Xs = 3*105 and 1.2 < p e /p 0 £1.4. Comparison 

with experimental data by Hakkinen et al [4l] makes it clear 
that, as far as the division of pressure, the friction along 
the length of the plate and the length of the separated region 
are concerned, better results were obtained than from the 
methods applied previously. Both investigations achieved 
approximately the same accuracies, with some advantage on 
MacCormack’ s [67] side. Because of the lack of computed results 
for higher Mach numbers and more intense shock waves, no 
estimate of accuracy can be made for those ranges. 

The MacCormack method has been recently perfected by 
MacCormack, Baldwin [ 69 ] (1975) and was adapted for use in tur- 
bulent flows. The grid form was also changed in that study 
so that Ay is changed through an exponential parameter which 
permits its outward expansion (similar to the Skoglund, Gay 
[94] method). 

In addition to these studies, two more are worth mention- 
ing, which deal with the numerical solution of complete 
Navier-Stokes equations for near (wake) problems. Allen, 

Cheng [2] (1979) solved for the field of flow in the wake of 
the rear step and Carter [13] (1972) solved the field of 
flow around the compression corner. Those two studies are 
based on the difference scheme of Brailovskaya [8] (1965) 

(which is also an explicit two-stage method). Allen, Cheng 
[2] made some changes in this scheme, which removed the 
Reynolds number from the conditions of CFL stability and 
caused acceleration of convergence particularly at low Reynolds 
numbers. In both of these studies the outer boundary is de- 
fined as very close to the end of the boundary layer, which 
is made possible by formulation of the boundary conditions 
for the external flow with uniform flow characteristics along 
the characteristic directions that extend from the field of 
flow to the outer boundary. Thus the influence of expansion 
and compression stream lines, that travel from the boundary 
layer to the outer boundary, can be considered. 

2. 1.3. 4. Summary and Comparison of Experimental and 
Theoretical Results 


Figure 5 presents a comparison of results from the princi- 
pal systems demonstrated with experimental results (the equation 
was taken from Murphy [74] and MacCormack' s [67] results were 
also plotted on it). 
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Summing up the results : /23 

a) The methods that employ the complete Navier-Stokes 
equations provide better solutions than those obtained from the 
boundary layer methods (at least within the range of flow con- 
dition under which they were obtained, i.e. M Q = 2; Re Xs = 

= 3*10^ and 1.2 < p /p < 1.4). 

G U 

b) For high Mach numbers and intense shock waves only 
results from b'oundary layer methods are available and those 
diverge the more from experimental results the larger the Mach 
numbers and the more intense the shock waves are. It appears 
that some more improvement is desirable for the solutions ob- 
tained with Navier-Stokes equation methods, as well as the 
extension of their use for a broader range of parameters Pr; 

Re Xs ; P e /P 0 ; M q . The influence of heat transfer (heating 

and cooling) and mass transfer (suction and injection) along 
the plate should also be investigated. 

2. 1.3. 5- Methods for Analysis of the Separation Phenomenon 

Since the phenomenon of boundary layer separation is tied 
directly to interaction and is a measure of its intensity, it' 
is important for practical considerations. We will (in addition 
to solutions of interaction that were demonstrated in previous 
paragraphs) here discuss several studies that deal with the 
analysis of separation against the background of various 
problems • 

2. 1.3 -5.1. Theoretical Analysis of the Separation 

Tani [ 98 ] (1953) and Page [77] (i 960 ) discussed the pre- 
diction of separation with formulations of approximation (based 
on boundary layer equations). Stewartson, Williams [97] (1 969) 
thoroughly analyzed the separation of flow during interaction 
between an incident shock wave and the boundary layer by means 
of equations (constructed from asymptotic developments) 
approximated for three different regions: the inner boundary 
layer, the outer boundary layer and the separated region, 
through adaptations of interregional connections. An estimate 
of the type of interference needed to create the separation 
and of the size of the separated region, as a function of the 
flow conditions, is gained from the results of the analysis. 

Additional development of methods for prediction of 
separation from the intensity of the pressure gradient was 
carried out by Gerhart [33] (1973) and Messiter et al. [71] 

(1971). 
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2 . 1 . 3 - 5 - 2 . Practical Prevention Methods 
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A detailed description of flow separation phenomena and 
their negative consequences on aircraft performance is given 
by Cooke, Brebner [17] (1959)- Pearcey [8l] (1961) and Wuest 
[105] (1959) describe separation phenomena caused by various 
types of interaction (particularly for an incident shock wave 
and in a compression corner) and discuss different ways for the 
control of the boundary layer to prevent separation. 

The principal ways suggested are: Suction at the surface 
of the plate, special geometric design of the surface of the 
plane and artificial transition to obtain turbulence of the 
boundary layer where the capability of resistance to the pres- 
sure gradient is still stronger than in the laminar condition. 

Gai [29] (1969) discussed the experimental results of 
the influence of suction on interaction between the shock wave 
and the turbulent boundary layer in the compression corner. 

Anderson et al. [3] (1969) describe a combined method of 
suction in the separated region and injection in the stagna- 
tion region, for improvement of laminar characteristics and 
reduction of friction on a rough surface with a high Reynolds 
number. 

Inger, Swean [48] (1975) analyze the influence of suction 
or injection in the deflection corners from the plate to the 
flow on a laminar boundary layer with a separated region (the 
analysis is performed through an approximated solution based 
on the Levy-Lees transformation). 


2.2. Finite Difference Methods of Solution of Navier-Stokes /25 
Equations 


This part discusses finite difference methods that were 
developed in recent years to solve Navier-Stokes equations 
for many and varied flow problems. 

2.2.1. Background References 

The principles for the solution of differential equations 
by means of finite differences were first discussed at length 
by Courant, Friedrichs, Levy [18] (1928) who classified the 
types of problems and the various methods of computation and 
prepared tests for the stability of the solution. Lax [57] 
(1967) surveyed their approach with reference to recent devel- 
opments on the subject of computation with finite differences. 
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Mathematical background, general numerical methods and 
discussion of the stability of the solution through application 
to practical problems are to be found in books by Forsythe, 
Wasov [27] (i 960 ) and Richtmayer, Morton [ 85 ] (1967). 

A detailed survey of numerical methods for all areas of 
fluid dynamics is given in Roache's [86] (1972) book. 

A general new formulation of various difference methods 
and the connections between them in the solution of all types 
of differential equations is given by Isenberg, Davis [49] 
(197^) where the principles of consistency, stability and 
convergence are also discussed. 

2.2.2. Finite Difference Methods 


The approach to the solution of Navier-Stokes equations 
through finite differences was thoroughly outlined in the 
investigations by Crocco [20] ( 1965 ) and Cheng [16] (1970), 
who suggested a way for an asymptotic solution of the Navier- 
Stokes equations that are functions of time, which aspires to 
a solution in the constant condition. The process of inte- 
gration over time is analogous to the iterative process for 
the solution of equations that are not functions of time. That 
way is justifiable since it is a more correct model of the 
physical process. Cheng [16] discussed in great detail the 
stability of the solution and the need for introduction of 
simulated viscosity by means of a cutoff error in the differ- 
ence equations for the sake of a stable final solution. 

The equations are formulated in the form of a group (by 
means of variables p, pu, pv, E) since that way the gradients 
of the variables are smaller and there is less chance for 
interferences with the numerical computation that might develop 
because of them. 

Of the existing methods we will discuss below principally 
those that deal with the solution of equations that are time 
dependent. In the studies by Lax, Wendroff [ 58 ] (1964), 
Brailovskaya [8] ( 1965 ), as well as Palumbo, Rubin [ 78 ] (1972), 
basic procedures that are very much alike are discussed for 
the solution of a series of equations in their [term unknown] 
form through an explicit two-stage method. These methods 
differ from each other only in their choices of technique to 
compute the finite differences at the grid points. The tech- 
nique chosen has significant influence on the cutoff error 
that serves as artificial viscosity and plays an important 
role in stabilizing the solution for the series of equations. 



Based on the Lax, Wendroff [58] method, many variations 
were developed for various cases. Rubin, Burstein [90] (1967) 
used this method in solving for inviscid fluids with shock 
waves. Thommen [100] (1966) investigated the stability of 
viscous flow on a plate and around a step with it. 

MacCormack [67] ( 1969 ) developed perfections and improve- 
ments to the basic scheme and used them to solve the interaction 
problem in MacCormack [71] (discussed in this research). 
Skoglund, Gay [94] (1969) also solved this problem by introduc- 
ing some changes into the scheme. (The last mentioned papers 
were discussed at length in Section 2. 1.3. 3-) 

Even the Brailovskaya method [8] was adopted for the solu- 
tion of various problems, the most important of which are those 
of Allen, Cheng [2] (1970) who solved the flow behind a step, 
and of Carter [13] (1972) who solved the flow around a compres- 
sion corner. (These papers were also described in detail in 
Section 2. 1.3. 3. ) 

The Brailovskaya [8] and MacCormack [67] schemes were 
adapted for use in this research and are formulated in detail 
in Part 4. 

The Palumbo, Rubin [78] method is similar in principle to 
those of Lax-Wendroff and Brailovskaya. Its distinction lies 
in the better arranged computation of differences, according 
to which the calculation of each point depends on a larger 
region of points (12 in number) around it. This method solved 
the problem of compressive flow before and after a step. 

The Harlow, Amsden [43] (1971) method is basically differ- 
ent from the previous ones. It describes an implicit scheme 
for two-dimensional or three-dimensonal flow in which the 
pressure is calculated from Poisson’s equation (which is 
obtained from the momentum and continuity equations). In this 
method the variables are solved as non-identical points in 
the grid meshes (p, p, T and u, v are calculated in two com- 
bined grids alternately by the Staggered Grid Method). 

The Scala, Gordon [91] (1967) method is built on an 
implicit-explicit scheme, using the exponential transformation 
of the coordinates. There are two schemes to this method: one 
is implicit and solves the unknowns in each field through a 
system of linear equations that are solved Iteratively; in the 
second, explicit scheme, the unknowns are computed separately 
for each point by means of repeated substitution. The two 
schemes are used alternately. 

The last three methods noted are relatively new and have, 
in the meantime, been tried out for very simple problems and 
only for occasional special cases. 



It seems proper to take note of some additional studies 
when comparing results of various numerical methods as to 
their accuracy and stability of solution. Brailovskaya et al. 

[9] (1970) and also Dorodnitsyn et al. [22] (1972) review 
developments of various methods, particularly for incompressible 
flow. Emery [24] (1968) and Tyler et al. [99] (1972) compare 
a number of methods for the solution of viscous flow problems 
and of a one-dimensional shock wave transition in inviscid flow. 

2.2.3. Computation Techniques for Flows Including Shock Waves 

Due to the physical (and mathematical) discontinuity of the 
shock wave strong inhibitions to a numerical solution arise, 
which necessitate the employment of varous numerical techniques 
for their prevention. The most accepted one is the use of 
artificial viscosity through the cutoff error of the scheme or 
by addition of more terms to the scheme. The artificial 
viscosity "spreads" the discontinuity of the flow over a wider 
region than the actual physical width (this method is known as 
shock capturing). In all of these methods the artificial 
viscosity is so constituted as to be proportional to the 
gradients in the field of flow. 

The following studies are concerned with quantitative and 
qualitative evaluations of artificial viscosity and with 
recommendations for its use in different cases. Hirt [44] 

(1968) investigated the stability of nonlinear equations based 
on testing the cutoff error, Davis, Mallinson [21] (1972) 
tested the influence of the cutoff error from preceding differ- 
ences in various flow problems, in computation grids of a 
different size. Roache [87] (1972) found a connection between 
the various influences of artificial viscosity on the solution 
of the implicit scheme in the transient and steady states. 

Goodrich et al. [37] (1972) developed a perfected method 
for adding artificial viscosity terms to the difference scheme, 
which is based on consideration of the gradients in both 
dimensions of each field of flow. (This method was used in 
the present research and is described in detail in Part 4.) 

Another way to overcome discontinuity was modelled by /28 

MacCormack, Warming [68] (1973). The shock wave is removed 
from the region of the finite difference computation by 
assigning double values (which are connected to each other 
through shock wave transition formulas) to grid points passed 
by the shock wave and by calculation of preceding or following 
differences in their vicinity so as to avoid calculation of 
derivatives for the shock wave path itself. This method is 
called shock splitting and it was also tried (with partial 
success only) in this research. 



2.2.4. Finite Differences for Special Cases Related to 
Research Subjects 

2. 2. 4.1. Boundary Layer Solutions 

Since the methods for the solution of boundary layer equations 
with finite differences are important for understanding the solu- 
tion of Navier-Stokes equations, we will describe their sources 
briefly. A detailed survey of existing difference methods is 
given by Paskonov, Chudov [80] (1968) and Blottner [7] (1970). 
Blottner [7] and also ‘Waiter, Leblanc [102] (1971) and Jaffe, 

Smith [50] (1972) developed methods for the computation of 
boundary layers that also contain equations of the system’s 
chemical composition. Jaffe, Smith and Lubard, Shetz [ 65 ] 

(1968) also discuss the influence of changes in the boundary 
conditions (line injections, suction, heating and cooling) on 
the results. The difference method by Fliigge-Lotz, Er-Young 
[ 26 ] (i 960 ) is particularly suitable for help in the solution 
of pressure interaction problems. Werle, Bertke [ 10 4 ] (1972) 
occupy themselves with computation of the separated boundary 
layer by assumption of a velocity profile and its introduction 
into the computational model. 

2. 2. 4. 2. Approximate Solutions of Navier-Stokes Equations 

These solutions, even though they are not applicable for 
the most general problem, are of significance for cases where 
the entire field of flow, or part of it, is to be described 
approximately and in simple fashion. 

Some solutions of Navier-Stokes equations through finite 
differences for viscous and noncompressive flow are given by 
Loer [64] (1971); Lascaux, Raviart [ 56 ] (1970); Welch et al. 

[103] (1965) and also by Ghia, Davis [34] (1974). The non- 
compressibility permits transformation of the variables u, v, 
p to 1 ) 1 , in through release of pressure, which enhances the 
final solution considerably. Gosman et al.’s [39] (1969) book 
describes the method of general solution for noncompressive 
flow with a change in heat, mass and chemical composition. 

Jamet [51] (1970) discusses the difference method for compres- 
sion equations that are inviscid and time dependent. The study 
by Gottlieb, Gustaffson [40] (1974) assumes an approximate 
algebraic expression for the energy equation by assumption of 
a fixed overall temperature and solves the hyperbolic flow 
through finite differences with artificial viscosity. 


2.3. Summary of the Main Changes and Improvements in the /29 

Present Study 


Below we will summarize the main changes and improvements 
to be found in this study, with reference to the most advanced 



existing studies; we will also describe the principles of the 
experimental changes that were made, but not included, in the 
final computation methods for various reasons. The studies to 
which this part refers are those by MacCormack [ 67 ] and 
Skoglund, Gay [94] (which solve an interaction problem similar 
to the subject of this study), as well as by Carter [13] 

(solution of interaction around a compression corner through 
characteristics in the outer boundary). 

2.3-1. The Computational Scheme 

The basic difference scheme was prepared according to two 
methods in parallel: the Brailovskaya [8] method (according to 
the Carter [13] formulation) and the MacCormack [ 67 ] method. 

A number of changes and options were tried out with the basic 
scheme, to be used in accordance with specific purposes: 

a) Change of grid structure in the x direction (not included 
in the final scheme) : 

In addition to the regular, uniform grid (where Ax/Ay = 16) 
there is a possibility to use a nonuniform grid, dense in the 
middle and spread out toward the corners (where 1.6 < Ax/Ay < 50 
and the ratio between the lengths of two adjacent meshes is 1.1), 
similar to the Skoglund, Gay [94] method, which is detailed in 
Paragraph 3 -4. 3- 2 and in Appendices H and I). This option was 
not included in the final scheme since it caused large trunca- 
tion errors and strong oscillations in the solution for the 
vicinity of the shock wave region (as detailed in Paragraph 
4.1.1). 

b ) Differences of the fourth order in the x direction under - 
neath the shock wave (not included in the final scheme): 

This possibility helped at first to increase accuracy of 
computation in the regions with sharp gradients in the x direc- 
tion (in a grid where Ax >> Ay), when the differences in the 

x direction are computed to an accuracy of O(Ax^), as against 

the regular computation with an accuracy of 0(Ax2). Details of 
this method are given in Paragraph 3- 4.3.1 and in Appendix H. 

This option was also not included in the final scheme since 

strong oscillations, whose origin is in the contact area between 
the regions of second order computation and fourth order compu- 
tation, were noted in the computation. (Additional details are 
found in Paragraph 4.1.2.) 

c) Addition of artificial viscosity (included as option 
in the final scheme): 

The addition of terms for the artificial viscosity aids in 
stabilizing the solution, particularly for intense interactions. 
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Its employment helps, particularly during the first stages of 
computation when the results are still far from the final solu- 
tion and subject to great changes. Expressions for artificial 
viscosity are based on those by Skoglund, Gay [94] with slight 
changes, as detailed in Paragraph 3.4.4. The use of this option /30 

aids in the solution but requires preliminary computation for 
optimization of several numerical coefficients (details in 
Paragraph 4.1.3)- 

d) Testing stability and . convergence (included as option in 
the final scheme): 

The possibility is offered to calculate the relative change 
of the variables u, p, T and the remainders of the four differ- 
ential equations for each point of the computational grid, at 
each step of the iteration. The average maximum values of these 
changes are a measure of the stability of the computation and 
the convergence of the solution (as defined in detail in Para- 
graph 3- 5- 3-1). These tests are not performed automatically 
during the computation; they are rather offered as optional 
choice (because they increase the memory stored and the computa- 
tion time by nearly a factor of two). 

2.3.2. Boundary Conditions 

a) Removal of the shock wave from the computational field 
(not included in the final scheme): 

In addition to the usual difference computation over the 
whole field of flow with the shock wave, the possibility was 
provided to remove the shock wave from the field so that the 
flow conditions in front of it and behind it would serve as 
boundary conditions, as described in Paragraph 3- 3- 1-5- This 
method (which, in the past, has not been tried for this type 
of problem) is called shock splitting, or the method of "dis- 
continuous" computation of the shock wave path. The use of 
this method does not appear suitable at this stage since 
significant oscillations were generated along the graded border 
of the shock wave (as explained in detail in Paragraph 4.1.4). 

b ) Boundary condition of the external flow characteristics 
(exists in the final computation scheme): 

This boundary condition is defined similar to Carter's [13] 
formulation and its main purpose is consideration of the influ- 
ence of streamlines from the reflected shock wave on the outer 
boundary. This permits the cutting down of the computational 
field to the region close to the termination of the boundary 
layer (in contrast to previous studies in which the outer border 
was extended so that the reflected streamlines would not exit 
through it but rather in the downstream direction, a procedure 
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requiring a computational grid that is much too big). The method 
is explained in detail in Paragraph 3. 3. 1.4. 

c) Boundary condition for density along the plate (exists 
in the final computation scheme): 

This condition is computed (due to considerations for the 
stability of the numerical computation) on the assumption that 
at close range in front of the plate, all along it (in the 
separated region as well) there prevails the condition 

U.\ -o ■ t>P| a i , 

•jjvj \ w i as well as the adiabatic condition (as ex_ 

plained in Paragraph 3.3.I.2.). In some of the previous studies 

this boundary condition was computed from the solution to one 

of the equations for flow next to the plate, or through extrapola- /31 

tion of the density from the field to the plate. Such procedures 

are essentially more accurate but are sometimes cause for 

instability in the numerical calculation; the above-mentioned 

method is therefore preferable. 

d) Polhausen profile at the entrance cross section (exists 
In the final scheme): 

Boundary conditions at the entrance cross section are com- 
puted according to the Polhausen method for profiles of flow 
in supersonic boundary layers (Paragraph 3- 3. 1.1). 

2.3.3. Parametrical Analysis of Interaction Characteristics 

The computations described in this paragraph were carried 
out according to the final scheme in which the grid is uniform 
in the x and y directions; second order differences are applied 
for the entire computation; the shock wave is inside the field 
("continuous" computation); artificial viscosity has been intro- 
duced and the boundary conditions b), c) and d) from Paragraph 
2.3.2 apply. 

a) Speed of convergence 

In the above-mentioned computation method the solution 
converges, after only a few hundred iterations, to results with 
a relatively large truncation error; yet their accuracy is 
only slightly inferior when compared to previous methods of 
computation (which required considerably longer computation 
time). This -way the computation is serviceable and practical 
and, indeed, 100 different computer runs were completed for 
testing the parameters of the results, each run taking up only 
3-5 minutes of CPU time. 
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b ) Parametric analysis of the results 




While in previous studies analyses of results were carried 
out only in a few cases (because computation time was relatively 
long) a detailed analysis of the interaction was performed in 
this study, according to the principal parameters that influence 
it. The influence of the Mach number in the range 2 < M Q < 4.5, 

the Reynolds number 10^ < Re Xs < 10^, the shock wave intensity 

1.2 < P e /p Q < 3.2, the shock wave distance from the leading edge 

3 cm < x < 9 cm and of the Prandtl number 0.72 < Pr < 1.0, was 

tested. Also tested was the influence of heat transfer (a non- 
adiabatic wall) and mass transfer (suction and injection) along 
the plate. 

c ) Detailed analysis of the influence of suction on the 
separated region 

A detailed analysis of suction influence, in combination 

with various flow conditions, on the interaction in general and 

on the separated region in particular was carried out. The 

required amount of suction, for each combination with Re x ; 

s 

P e /p 0 i M 0 to avoid boundary layer separation, was determined. 

These results have much practical application since the prevention 
of separation improves the performance of aircraft components — 
particularly of control surfaces and jet engine inlets — around 
which this interaction takes place. 


3. Computation System /32 

3.1. Basic Equations 

3.1.1. Dimensional Flow Equations 

The basic flow equations are the Navier-Stokes momentum equa- 
tions, the continuity equations, as well as the fluid state 
equation and the formula for the temperature dependence of visco- 
sity. The equations meet the following assumptions: 

a) The coefficient of bulk viscosity is neglected. 

b) The coefficient of viscosity is a function of temperature, 
according to the Sutherland formula. 

c) The state equations are for an ideal gas. 

Accordingly, the noninstantaneous equations are written for a 
viscous, compressible, heat-conducting fluid in a two-dimensional 
f Cartesian system (per Schlichting [92]) 
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Continuity equation . 




( 3-la) 


Momentum equation 
in the x direction 


(3-lb) 


«* + ( + pV£if = t(cr A Z)*±Ai*) 


Momentum equation 
in the y direction 


( 3-lc) 


* 




a/ . pV ^ + pV 1 ifxy) ■> -v(C^l) 


*3 


Energy equation 


(3-ld) 




Equation of state 


p% r /t 
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(3-2) 


Sutherland's formula for viscosity 


f V= c svr v 


(3-3) 


*This assumption is correct for a monatomic gas in which there 
is no internal degree of freedom for the molecule. This does 
not hold true for a multiatomic gas and the coefficient can 
reach the order of magnitude of the regular viscosity coefficient, 
particularly within a shock wave. Since there is no good resolu- 
tion within the wave in the present computation the volumetric 
viscosity can be neglected. Additional details on this subject 
are offered by Vincenti, Kruger [101] (1965). 
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when 


ftciir =287 .04 - 77 ^ 


sec'- a < ! 


S* =110. ^ K. 


Stresses In the momentum equation are defined by 

« • « * '*• 1 - i"“1 ■ jK(*S-iS) 


}?(*}£■&) | « 

I 

r \ ^ m * tjx */ 


See note at the end of the paragraph. 

The dissipation function In the energy equation is defined by 

* - r1-$Sf-«S i(8-?S)'J 


( 3 - 5 ) 


The terms for heat conductivity are given by 


G *’ ic 


( 3 - 6 ) 


ef « v* n: 


We define the function of total energy for the energy equation as 




( 3 - 7 ) 


Note: The stresses in the momentum . equation are based on Hooke’s /3^ 
Law for fluids, whose general formula is 
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SfT • (- ) a + fx" (As - A, - -) 


Stokes' assumption jh + -'o is set against this formula 

(its significance being the neglect of volumetric viscosity). 
That is how we obtain Eqs.ji';(3-4) . 


3.1.2. Definition of Non-dimensional Variables 

We will define the system of non-dimensional variables as 
follows : 



^ 4 t. ±1 « ii 

u. . 


v*: 

uj- 

f-'4 

P» P * T- T* 

j 3 . 

uV-/Cp 

> 

ii 

3k 

y- £e- 5 . s* 

H = 

= y[ T /V + 


*■* 


A 


V 


J 


( 3 - 8 ) 


We will also define the coefficients: 


Reynolds number 


Mach number 


Prandtl number 




M = 


r rt 4 - J 



p,.-- 


K* 



( 3 - 9 ) 
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When x* s - 

<>o - 


- reference length (distance from the leading edge to 
the shock wave continuation impinges on the plate); 

- condition of external flow at the entrance cross- 
section. 


We thus obtain the non-dimensional values for external flow 
at the entrance cross-section (reference values): 


* { 


(H) M . 1 


(3-10) 


J 3 * ~ ^ p' j - 4 


Since (V#) = 0 (see Paragraph 3* 3* 1.1), we obtain V 0 = 0 when 




3.2. System of Non-dimensional Equations 

A system of non-dimensional equations is obtained from these 
last expressions. The four differential equations will be written 
in group form for the group variables p, pu, pv, E. This form 
is handier for use in numerical computation (the gradients of 
these are smaller and prevent the tendency to diverge). 

The general form of the differential equation is 


. 3 F & __ ^ 


— — + — 
TiX 


(3-11) 


when 


P<-pu L 

jpuv 

U(E^P) 


«- x y 


* 


f* 

S'** 

P+fV 1 - 

V(£±f) 


' lA."t"xv Jr *'■ 


( 3-lla) 
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when (in non-dimensional formulation) 


/36 



the equation of the non-dimensional state is 


(3-llb) 


P = 



The formula for non-dimensional viscosity is 


(3-12) 



5 + 


* 3/ 

{r-0 Hg t, t /(ar-«) M^t 


4 * S( 

4 + Vr 


(3-13) 


for S = 110°K 



5* 

fr-l) H. v T 0 fr I 


since the formula for non-dimensional viscosity also depends on 
the temperature (dimensional) T* Q of the external flow. 

Four differential equations (3-11) and two algebraic equa- 
tions (3-12) and (3-13) represent the system of six equations 
in six unknown quantities, when the procedure for solution is 
in fact carried out through the differential equations (into 
which pressure p and viscosity y are inserted from the algebraic 
equations) from which solutions for the unknowns p, u, v, T 
are obtained. 


In addition to the solution for those unknowns, the Mach 
number and the flow function in the field of flow are also 
computed by means of the following non-dimensional connections: 


Mach number 



(3-14) 
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Flow -function 


(3-15) 






when a, b 'are random points in the field of flow (the 
function is computed according to the definition 



^ & 


3 


fv'-~ 


1st 

x 


flow 


The coefficients of friction and heat transfer along the 
length of the plate are also computed by means of dimensionless 
values , as follows : 
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Friction coefficient 

Heat transfer 
coefficient 


-i. 

fte- I » 




w 


fr-dHf.i'H. 


w 


( 3 - 16 ) 

(3-17) 


(Computation of heat transfer is carried out only for non-adiabatic 
boundary conditions.) Details about development of the last four 
non-dimensional equations are given in Appendix A. 


3.3. Initial and Boundary Conditions 

The above-defined equations constitute a system of partial 
differential, nonlinear equations of parabolic characteristics 
vs. time, which, in its permanent form, becomes a system with 
aliphatic characteristic. The problem is defined as an initial 
value-boundary value problem and requires definition of the 
boundary condition along the length of all boundaries of the 
field of computation and the initial condition inside each field 
of computation. For the four differential equations we require 
boundary conditions and initial conditions for the components 
of velocity v; u and for density p and temperature T. (The 
pressure p and viscosity y are computed, according to the 
boundary and initial conditions, by means of the algebraic 
equations of state and viscosity. ) 


3.3.1. Boundary Conditions 

The boundaries of the computational field include four 
sections, as defined in the following diagram. 
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AD Entrance cross section 
AB Plate 

BC Exit cross section 
CD External flow 

E 2 F 2 (as shown in Fig. 2). 


According to the compu- 
tational method where the 
shock wave is removed from 
the region considered by 
the computation (the "non- 
continuous" method), the 
external flow included two 
sections F^C and DE^ plus 

the boundary section up- 
stream from the shock wave, 
E -^2 and the boundary sec- 
tion downstream of it, 

E 2 F 2 (as shown in Fig. 2). 


3. 3. 1.1. Entry Cross Section Conditions 738 

With the assumption of and *^-1 ' =o in the vicinity 

of the entry cross section, a velocity profile is obtained through 
Polhausen’s method (which assumes the velocity profile as the 
progression range of the coordinate normal to the plate, 
which fulfils the conditions on the plate and at the edge of 
the boundary layer) in the following manner (all expressions 
non-dimensional) : 


W- . \ 

ti = * *1- 


U.O 

UoM 


[,-(*)'] j. 

T ‘ ’ a-ijHJ' J 


(3-18) 


(3-19) 


Under assumption of fixed pressure (since 3p/3y = 0 at the entry 
cross section) 

l - +7 

f* T/To 




( 3 - 20 ) 


The profiles of T; u; p are here given as a function of n> 
the connection between n and y is given as 




(3-21) 
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Component v of the velocity at the entry cross section may 
be neglected 


V = 0 (3-22) 

Detailed development of the boundary conditions at the entry cross- 
section are given in Appendix B. 


3. 3.1.2. Boundary Conditions Along the Plate 

From the nonslip conditions along the plate, we obtain 
for velocity components 


u = 0 (3-23) 

v = 0 (3-24) 

When mass transfer via the plate is made possible, then the v /39 

condition changes to: 

For injection v > 0 

(3-24a) 

For suction v < 0 

(In these cases the value of v is on the order of a few percent 
of u. ) 

The adiabatic condition along the plate dictates the boundary 
condition for temperature T 


(3-25) 

i 

When the plate is kept at a fixed temperature, the T condition 
changes to 

T = const. (3-25a) 

The boundary condition for density can be computed in several 
ways : 


(a) from one of the motion equations (continuity or momentum), 
which assumes a simpler form next to the plate; 

(b) through extrapolation from the third or fourth order 
of density in the field of flow, towards the plate; 

(c) from knowledge of pressure and temperature at the plate 
(through the state equation) when the temperature is given 
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per (3-25) or (3-25a) and pressure at the plate is computed 
on the assumption that next to the plate the relation 

lL*o \ 

■ ( 3 - 26 ) 


is valid. This condition (which prevails in the regular boundary 
layer over its entire width) is particularly well fulfilled in 
this interaction in the lower part of the boundary layer (which 
is close to the plate) and in the separated region, in spite of 
the interaction influences. The assumption is therefore justi- 
fied (see the note). 

Of the three ways available, the last one was chosen due to 
consideration of stability for the numerical computation (even 
though it is less accurate, basically, than the other two). 

The main reason was the tendency of solutions available through 
the other two methods (particularly the first one) to oscilla- 
tions, since the local computation of the derivatives close 
to the plate is not accurate and distorts the numerical values. 

Note: During previous tests in this investigation (before /40 
formula ( 3 - 26 ) was introduced for computation of the boundary 
conditions) and also in other computational studies, it was 
z' found that the assumption of 3p/3y = 0 is valid as a very good 

approximation for the entire length of the plate; its use is 
thus entirely justified. 


3. 3.1. 3. Exit Cross Section Boundary Conditions 

The location of the exit cross section is so determined that 
no further significant changes in flow characteristics can occur 
in the x direction. It can therefore be assumed with good 
approximation that 


bu _ 'dv _ 

7>fc" / dX~'&X 'OX. 


(3-27) 


That is a bit misleading, because there are very small 
gradients in the x direction (because of the renewed development 
of the boundary layer downstream of the Interaction), but due 
to the hyperbolic characteristics of the boundary layer these 
distortions in the flow field (which are relatively small) do 
not have any influence on the interaction region; so the 
distortion, if it exists at all, is only a local one. 
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3. 3-1. 4- External Flow Boundary Conditions 

The location of this boundary is so determined as to be in 
the inviscid region, as close as possible to the plate but 
sufficiently far from the region where changes in the interaction 
are taking place. The order of magnitude of this distance is 
( 2—3 ) <5 o (a computational assessment of this magnitude is given 
in Appendix D). 

In contrast to previous studies (like Skoglund, Gay [94] 
and MacCormack [67]) where this boundary was chosen to be 
sufficiently distant so that its boundary conditions were 
uniform and fixed (and reflected waves- went out through the exit 
cross section), in this investigation the boundary conditions 
are influenced, because of the proximity of the external 
boundary to the region of interaction, by the form the develop- 
ment of the boundary layer assumes and by the compression and 
expansion streamlines that emanate from it. These conditions 
are thus not of uniform nature and they must be computed 
point by point from the actual conditions in the flow field. 


Uniform flow characteristics are therefore used for flow 
outside the boundary layer along the directions of the charac- 
teristics, so that the external boundary condition is updated 
during the computation process in accordance with changes of 
the flow characteristics in the field. 


Basically, this is done in 
the following stages: 



(a) A reference line is 
fixed from which characteris- 
tics* fan out to the external 
boundary; this line is close 
to the external boundary and 
outside the boundary layer; 


(b) For each point A on 
the reference line the direc- 
tion of the local flow A0 and 

of the local Mach angle y is computed; their sum 8 gives the 
direction of this characteristic; 


r ^ ( 



( 3 - 28 ) 


/4l 


A 

(c) Point B is fixed at the external boundary (where the 
flow characteristics are identical to those of A) by extending 


*Two characteristics emerge from each point, as shown in the 
diagram, but for the purpose of showing the influence of the 
streamlines on the external boundary, the upper characteristic 
is necessary (which reaches that boundary). 
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a straight line from A to the external boundary at angle g. An 
approximately straight characteristic* is assumed, which is 
justified since Ah is very small. This computation is done 
along the length of the external boundary, except for the point 
where the shock wave enters the field. 



(d) Computation for both 
sides of the shock wave at the 
external boundary is carried 
out by means of Rankine’s 
and Hugoniot's shock wave 
formulas. The shock wave 
is defined as pressure jump 


r~- 


h Pe * 

which is compelled at the external boundary between the two 
reference points E and F. The flow characteristics at point F 
are computed from those of E (where the values are fixed by 
the characteristics scheme explained in the foregoing paragraphs). 

Shock wave transition formulas for stage (d): 

= (u e ^c r+v-t^r) c ^ — (f -v (ut^ricr- i/ g £ c.er^cri ( 3 _ 29 ) 

' (>-0 
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(U 


CT4- / u ^(3-30) 


T* = T £ 


( 3 - 31 ) 




iu_ v cT 


j. ( r-*) 

The pressure given is provided by the compulsory pressure jump 


fV -- • Sk 

Density is provided by the state equation 


( 3 - 32 ) 


*The characteristics should actually be curved because straight 
characteristics are present in non-rotational flow only; this 
non-rotationality is not preserved in this instance because of 
the intersection between the shock wave and the streamlines, 
which causes curvature of the wave. This phenomenon is neglected 
when the computational grid is very fine and so it is permitted 
to assume the characteristics of this region to be straight 
lines . 
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(3-33) 




The angle of the shock wave a is given by the relation 


(3-3*0 


Details of the development of these equations are given in 
Appendix P. 

It must be pointed out that since the flow characteristics 
upstream of the shock wave at E change as functions of flow char- 
acteristics in the field, the characteristics of P will also be 
dependent on these changes so that here is an expression of the 
influence of interaction on the incident shock wave itself (while 
previous studies assumed fixed values for both sides of the 
shock wave for the entire process of computation). 

Note: In general the pressure ratio p e /p Q serves as charac- 
teristic parameter for the interaction, and a ratio Pp/P E must 

be found that will match it. This procedure is iterative and 
explained fundamentally in Appendix C. 


_l 
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3.3«1«5» Shock Boundary Conditions for the Shock Fitting /4 3 
( "Noncontinuous" ) Method 

This method is intended to take the shock wave out of 
the field of computation, with the purpose of attenuating the 
disturbances of stability that arise because of the strong 
gradients along the shock path. 

In this method the boundary of the field of computation is 
assumed to be on both sides of the shock wave along its entire 

to the sonic line). 

The process of computa- 
tion is carried out as 
follows : 

(a) From the shock 
wave angle (Eq. 3-3*0 > 
we find the initial 
location of the shock 
wave in the field, and 
from it we get the 
boundaries e = E-^Eg, 


length (for the external flow boundary 
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upstream of the shock wave and f = F^F 2 , downstream of it 
(assuming- a straight shock wave); 

(b) Variable values along e are determined from the upstream 
conditions while values along f are determined from the values in 
e, by means of formulas for the shock wave transition given in 
Paragraph 3* 3 .1.4; 

(c) After updating the variable values in the -entire flow 
field (through a difference scheme) the' local gradient of the 
shock wave is brought up to date through 



and according to it the location of the shock wave is updated 


X*( 




+3 W 3 ) 3 h 


(3-36) 


The computation is continued per paragraphs (b) and (c). 


3.3-2. Initial Conditions 


/iMi 


Prom the purely mathematical point of view no importance 
attaches to the way initial conditions are chosen because, if 
there is a single value solution (which is assumed) it must be 
an expression that is not a function of the initial condition. 
However, when we talk of a numerical solution for an iterative 
process, which represents integration over time, it must be 
taken into account that the closer initial conditions come to 
the final solution, the faster convergence will occur and 
reduce computation time. 

Initial conditions are chosen so as not to conflict with 
two main limitations: 

(a) It is undesirable for initial conditions to be too far 
from the solution since oscillations in the values of the varia- 
bles can develop during the computation procedure (particularly 
during the early stages), which can not always be restrained 
(see also Roache [86]). 

(b) Initial conditions that are too close to the solution 
impair the general nature of the computational scheme and its 
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ability to solve various problems where the solution is not 
always known beforehand. , 

Initial conditions for computation have therefore been 
assumed as combination of the following data (as in Skoglund, 
Gay [94]) : 



(a) The profile of the boundary layer was assumed to be of 
slightly larger shape up to the shock wave and smaller again 
afterwards in the downstream direction, with a small separated 
region existing at the end of the incident shock wave; 

(b) The shock wave itself spread to the width of two meshes 
outside the boundary layer and to the width of 4-10 meshes in 
the subsonic part of the boundary layer. The flow characteris- 
tics downstream of the shock wave (in the supersonic region) 
were assumed to be uniform and were computed from formulas for 
shock wave transition. Definition of these conditions is simple 
and general and can be applied in identical manner for the 
entire combination of flow conditions. That way we get an 
initial guess at the solution, which allows speeding up of 

the convergence. 


3.4. Method of Finite Difference Solutions /4g 

3.4.1. General 

As mentioned before, two difference schemes were prepared 
simultaneously; that of Brailovskaya and that of MacCormack 
(though most of the computations were later carried out with 
the Brailovskaya scheme, since it offered savings in computer 
time and simplicity of operation, as will be explained in 
Paragraph 4.1.1). 

The two methods share a similar truncation error: 
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( 3 - 37 ) 


£r (mac.) = o( Ai v , ^ v ) 

£ T ( git A.) ” ° )’ 

V. 


The difference in dependence on At has almost no influence 
since the gradients vs. time are very small in the asymptotic solu- 
tion (as was also shown by Carter [131 in his comparison of accura- 
cies in various schemes). 


Those two methods have already been employed successfully in 
the solution for flows with mixed regions (parabolic, hyperbolic 
and aliphatic ones), as explained in paragraph 2. 1.3*3). 


3.4.2 Definition of the Region of Cbmputation 


Based on consideration for achieving the maximum possible 
accuracy and the minimum size of computational field required that 
would permit inclusion of all interaction phenomena (as is done in 
the sources listed in paragraph 2.1.3) the dimensions of the field 
of computation were determined as shown below: 


Length of field 
Width of field 


S'* 

i y c i 


r 


( 3 - 38 ) 


The considerations for this choice of dimensions are detailed 
in appendix D 1 . 

Location of the shock wave was always fixed so that the 
intersection between its continuation and the plate would be in 
the center of the region of computation. 

3.4.3 Computation Grid Dimensions /46 

3. 4. 3.1 General 


Difference computations were prepared simultaneously for 
two grids: one of them had a uniform mesh size while in the other 
the longitudinal dimensions of the meshes increased directly with 
their distance from the center. 
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The grid mesh dimensions were also set (similar to the field 
dimensions) as fixed multiples of i 0 and not as invariable longi- 
tudinal dimensions since the ratio between them and the boundary 
layer thickness is an important parameter which determines the 
accuracy of results. 


Computation of the number of points in the grid through 
these two methods is presented in appendix E’. Considerations 
for deciding the mesh dimensions are explained in appendix D’ . 

The various difference schemes in all methods, as well as an analysis 
of result accuracies, are detailed in appendices G* and H*. 

3. A. 3- 2 Uniform Computation Grid 
Dimensions for the grid meshes were fixed at 

1 

I (3-39) 

resulting in a ratio of Ax/Ay = 16. 

A number of difference schemes were prepared for this grid; 
two of them offer the accustomed accuracy (of the second order) 
with the Brailovskaya and MacCormack methods, one of them offers 
improved accuracy (of the fourth order in the x-dimension) with 
the Brailovskaya method, which is designated to increase the 
accuracy of computation in the regions with strong gradients in 
the x-direction. All of the schemes are explained at length in 
appendices G* and H*. 


3. 3. 3 Non-Uniform Grid 


Grid dimensions were fixed as 

^ i v. " o , "2, (To 

resulting in ratios 

. a V, _ ^ r 

= 1.6 If 


6 cfo 


“3 




(3-^0) 


Ax increased gradually towards the boundaries of the field 
according to the formula for geometrical progression, with a 
growth ratio of q = 1.1. 

3.4.4 Finite Difference Scheme 747 


Based on the general formulation of the differential equa- 
tions (3-11) 5 two difference schemes were written for the general 
case : 
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Brailovskaya scheme 
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with the truncation error .A^j ) 


MacCormack Scheme 
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For this scheme a division of the difference computation 
exists for the x and y directions. Detailed formulation of the 
two schemes and of the various stages of computation in them is 
given in appendix H' . 

The general procedure of computation is as follows: 

a) From the last values g£, (V| ; and 

are computed by means of definitions in formulas (3-11). 

. t-H 

b) \N vu. M is computed from the first stage of the differenc 

scheme and from it the expressions : , Q. and £ +1 

^ vi vu vi 

are calculated (but only for the MacCormack scheme). ' 

Vs/ 

c) >« ( vi is computed from the second stage of the difference 
scheme ana through it the final computation of the variables is 
carried out. 


The explicit definition of the solved variables in the equa- 
tion system is as follows: 





The variables determined 
by the differential 
equations 
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The variables determined 
by the algebraic 
equations 




f((r-0»CT 
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Artificial viscosity must be added to the difference schemes 
(3-41) and (3-42) by adding the expression C mjn to S mjn . (In fact 
this is only done for scheme (3-41)). This expression * is introduced 
into the difference scheme as an option, for the sake of stabilizing 
the solution during the early stages of iteration when the results 
are still far from convergence to the end solution. Sometimes 
this allows better control of the strong gradients, which makes the 
danger of divergence more remote. The definition of C m ,n is similar 
to that formulated by Skoglund, Gay [94] and Goodrich et al. [37], but 
it is somewhat simpler. 


C-K wstA 4. ^ V ** •' 
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when CCx and CCy are uniform numerical coefficients at all ' /49 
grid points (as against other methods, where the coefficients 
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change as function of the local Courant number). At any rate, 
these coefficients are chosen for providing optimum influence 
on the stability of the numerical solution. 

3.4.5 Boundary Conditions in Finite Differences 

The formulation of boundary conditions for finite differences 
in accordance with the definitions in section 3-3- is detailed 
at length in appendix J. 

3.4.6 Additional Computations in Finite Differences 

Computations of flow functions and Mach numbers in each field 
of flow, as well as of coefficients of friction and of heat transfer 
along the plate, are carried out with finite differences according 
to the definitions of section 3-3. Details of the computation are 
given in appendix J. 

3.4.7 Computer Program Procedure 

Fig. 4 shows a general flow diagram of the computer program 
as prepared for an IBM 168 computer. The region surrounded by the 
dash line includes the computation of differential equations and 
boundary conditions with finite differences. In the corners of the 
rectangles the names of the subroutines have been indicated. 

The diagram explains the basic computation for the Brailovskaya 
scheme, as explained in section 4.1. 2.8. 

All subroutines of the computer program, as well as the 
manner of employing the program in all its options, are explained 
in appendix Q. Also, a list of variables in the program with 
FORTRAN notations and details about their location in the various 
subroutines, is provided. 

3 . 5 Consistency, Stability and Convergence of the Computation 750 
Method 

3.5.1 Consistency 

The difference schemes defined in the previous paragraph 
contain truncation errors (which were computed in the reference 
sources in which the schemes were developed) of the following 
form 


Ar (llr-. i 1 ovskoya ) ~ o 
£7 OlacCormack ) ~ 


(3-45) 
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The truncation errors fulfill the consistency condition in 
principle 

(Hr) 0 

(3-^6) 

o 

The existence of this condition is given by the fact that 
the grid is fine enough so that all changes in the flow field 
can be identified in it. An extremely fine grid is, in fact, ' 
required for it, which is not practical. For that reason an 
optimum grid size is chosen, between desired accuracy of compu- 
tation and acceptable time for the computation and memory size. 

This interferes somewhat with maintenance of the consistency 
condition (as will be seen further on in paragraph 4.1.6). Di- 
gression from consistency is expressed by the final size of the 
grid meshes, causing truncation errors that are not negligible, 
particularly in regions with strong gradients (where Ay, Ax, At 
of the truncation are multiplied by high values for derivatives 
of the variables). This is expressed, among other things, by 
significant differences between the solutions for different schemes 
(for instance IVEacCormack [67] and Skoglund, Gay [94]), which 
appear to be consistency schemes on the fact of it, because of 
their truncation error. This phenomenon also stands out in dif- 
ferences between results received from the same scheme, when the 
grid is made fines (as in MacCormack [67]). 

3.5.2 Stability Criterion 

There is no possibility for demonstrating a criterion of 
stability that is both compulsary and sufficient for the difference 
scheme of the system of complete differential equations, because 
of their complexity and non-linearity. It is therefore customary 
to develop approximate differential equations and to employ them 
for the well known analysis by von Neumann, with its significance 
in testing whether the difference scheme of these equations 
restrains oscillations in the solution or increases them. The 
criterion received from this analysis gives the stability condi- 
tion for small disturbances only; this is only a compulsary 
condition (and not sufficient) in the general case since the non- 
linear expressions generally try to increase disturbances in the 
solutions and not to diminish them. The analysis relates only 
to testing reaction to periodic disturbances and also does not 
include a test of the boundary condition influences. The 
stability criterion is therefore only an approximate expression 
and in practical computations one needs to be more rigorous. 

Carter [13] prepared such a stability analysis for the 
Brailovskaya scheme along the following stages: 



a) Linearization of complete equations in their non-vectorial 
form for two cases - nonviscid flow (for which mixed derivatives of 
viscosity and the dissipation term were neglected) and viscous flow 
(for which the convection terms were neglected). 

b) Writing the difference scheme for the approximate equations. 

c) Placing of a single component of the Fourier series into 
the solution. 


d) Confirmation of the amplification matrix, which expresses 
the relation between values for variables in two adjacent iterations. 

e) Computation of the Eigenvalue of the amplification matrix 
and establishment of the criterion that would validate the 

von Neumann condition (which demands that this value be equal to 
or less than 1). Such a criterion is in effect a maximum value 
for At, because the Eigenvalue of the amplification matrix is a 
direct function of At. 


For nonviscid flow the condition 






- 


v3-47) 


is received (also known as the CFL condition after Courant, 
Friedrichs and Levy) while for viscous flow* 


St 4 


Pr' 




.At, 


(3-48) 
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is presented. In an analysis by MacCormack [67] of the 
scheme he developed, condition (3-47) is also found. 

MacCormack [ 67 ] and also Skoglund, Gay [94], who solved the 
system of complete Navier-Stokes equations did not at all refer 
to the condition At z of the viscous flow in spite of the fact 
that the following relation must be fulfilled 


A't - £ At a j (3-49) 


Z52 


The justification for it is found in. Carter's [13] analysis, 
who found that in such problems almost always Atq<At 2 so that in 
practice only the existence of At-, is required, in other words the 
CFL condition (3-47). 

Allen, Cheng [2] (who also solved the Navier-Stokes equations 
for flow in the wake of the final step) arrived at the conclusion 
that only the CFL condition must be validated by solving the system 
of equations for fixed viscosity, which prevented the existence of 
the condition for At 2 * 


* In other references expressions are found that are slightly 
different from (3-48) and they are given in detail in appendix K'. 
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Consequently, In this investigation the stability condition 
will from now on refer only to the CFL condition. Since computa- 
tion of the stability criterion is only carried out as approxima- 
tion this investigation will have to establish a more stringent 
condition, which is 

A-t = oC-.^ (3 _ 50) 

where a = 0.5 This value is in practice an optimum one since 
higher values brought the difference scheme to the threshold of 
instability, while lower values slowed the pace of convergence 
considerably. 

3.5.3 Test of Stability and Convergence of Results 
3. 5- 3.1 Variation Coefficients of the Results 


As measure of the convergence and stability of the numerical 
solution we have, in this study, employed several coefficients that 
are computed from the chronological results of the computation 
(in each area of iteration). 


a) 


Variation coefficients of the variables 
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b) Coefficients of the differential equation remainders 
(truncation errors) 
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when there are four systems for H for 

f V'; fly-* ^ ] 

(detailed development of the remainder coefficients is given in 
appendix L’). 

A maximum value and a mean value (r.m.s.) are computed for 
these coefficients in each field of flow. From the results it 
is possible to examine the degree of convergence by investigating 
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the development of the coefficients (their most frequent values 
or their averages) in the course of the iterations during compu- 
tation. Likewise, the division of coefficient values in the 
field may be checked and from it may be learned where the areas 
most sensitive to oscillations in the solution are. 

It is clear that to reach convergence these coefficients 
must decrease steadily down to the order of the truncation error 
for the difference scheme (or to the order of the computer's 
rounding error, if that is larger than that of the truncation 
error). 

Stability of the computation must be expressed in the monot- 
onous decrease of the coefficients' values as the iteration 
progresses, as well as those of their moderate and continuous 
gradients in the field of computation. 

3.5*3. 2 Estimate of the Order of Magnitude of the Compu- 
tational Error 

In the constant state* the Brailovskaya and MacCormack 
schemes have a truncation error of 0(AX 2 ; Ay 2 ) but it is very 
different to evaluate its exact numerical value since it is 
multiplied by various derivatives and those depend a lot on the 
nature of the field of flow and of the boundary conditions. 

According to the computation in appendix M' we get At = 0 ( 10— 3 ) ; 

Ay = 0(10- 3 ); Ax = 0(10- 2 ); and from them we learn that the fixed 
error of the variables (like for instance u, p, T) must be 
0(10- 6 ) (without consideration of the values for the derivatives 
with which Ax, Ay and At are multiplied and ignoring the boundary 
areas where the truncation error itself is considerably larger). 

According to appendix N' a computational estimate was carried 754 
out for the magnitude of the truncation error in a regular grid 
of 76 x 25 points, in comparison to computation for a grid twice 
as fine in each direction and containing 147 x 48 points. Accor- 
ding to this evaluation it was found that the error for the variables 
was on the order of 0(10- 3 ) . 

A check of the actual results of the computation reveals 
that most of the errors in the variables were in the order of 
0(10- 3 ) (per Pig. 10) and even the averages were around 0(10- 4 ) 

(per Fig. 11). 


*When the complete solution has oscillatory components there is an 
additional contribution to the truncation error, multiplied by At 
in the Brailovskaya scheme and by At 2 for MacCormack 's scheme. 



In the light of all this a final evaluation can be made 
for the computational accuracy of variables, which is in the 
range 0(10- 3 ) - 0(10- 4 ) . 


Note: The rounding error of the computer is in the seventh 

digit, which means the addition of an error of 0(10- 6 ) to the 
other errors of computation (because the definitive value of 
the non-dimensional variables is on the order of 0(1)). It is 
dlear that this error has no influence at all, as has even been 
proven in one computation for comparison with double the rounding 
error (of 13 digits), as explained in Fig. 14. 

3. 5. 3- 3 Comparison of Orders of Accuracies for Results 
with Different Methods 


In appendix 0* a table is shown comparing various methods 
for obtaining results to problems of interaction between shock 
waves and boundary layers. The table shows that the computation 
accuracy in the present study falls short by 1 -2 orders of 
magnitude of the published accuracies for other methods. This 
difference also exists between the original MacCormack method 
(in MacCormack [67]) and the same method employed in this study 
for special boundary conditions. This fact hints at the close 
possibility that the boundary conditions caused a larger compu- 
tation error (apparently because in the region of the incident 
shock wave, where the boundary conditions of the characteristics 
are close to those of the shock wave transition path, a signi- 
ficant local distortion is generated in the fulfillment of the 
flow equations, which also influences a part of the computational 
field that is not negligible). In spite of it this computation 
error does not affect the final accuracy of the solution very 
much, when compared to experimental results and previous compu- 
tational results (as will be made clear further on in part 4). 

On balance, the differences between results from all compu- 
tation methods (MacCormack [67] Skoglund, Gay [94], and the 
present research) and experimental results (Hakkinen et al [41]) 
are on the order of several magnitudes greater than the adver- /55 

tised accuracy of the methods of computation (particularly when 
close to the shock wave and in the separated region, which 
are two of the most important areas). There is therefore no 
great practical significance to the "lower" relative accuracy 
of computation in this investigation. 

Since a larger computation error requires less computation 
time for convergence, solutions were obtained relatively more quickly 
for this study than for previous ones (an additional factor that 
speeded the convergence was the choice of suitable initial 
conditions) and their accuracy, when related to experimental 
results, fell only slightly short of MacCormack’ s [67] work. 

(See section 4.2. 1.1). 
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This outcome shows, amont other things, that it is worth- 
while to accept larger computational errors as a tradeoff against 
shortened computer time, for the solution of problems of this 
nature. 



4. Results 


756 


4.1. Computation Method 


Prior to analyzing the results of computation with the 
difference scheme there is need to determine the structure of 
the final computation within which it will be executed (among 
the various suggested methods and options for computation). 
Considerations of accuracy of the results and simplicity of 
execution will be the determing factor. 

The equations that will be presented in the following para- 
graph all deal with the case in which 

R = 3 x 10 5 ; p /p = 1.4; M - 2 (below the reference 
e x g condition) 

This case was chosen because a separated flow already exists in 
it and there are various experimental results available for it, 
as well as theoretical ones, for comparison. The computer runs 
for this part are illustrated in Table 1. 

At first the basic schemes for the Brailovskaya and 
MacCormack methods will be compared, then the options and various 
changes made in the basic scheme will be investigated. 

4.1.1. Comparison of the Brailovskaya and .MacCormack Schemes 

The computation schemes were prepared according to the two 
methods by Brailovskaya [8] and MacCormack, as indicated in 
section 3*2 (except for the difference methods themselves, all 
other characteristic components were introduced into the compu- 
tation scheme, such as identical initial conditions and boundary 
conditions for both schemes). These schemes were applied in 
runs 001 and 051 equally and their results are compared in Figs. 

9 and 13 . 

It turns out that in the results from both schemes are very 
much alike, with a slight advantage to MacCormack’ s scheme 
(where the separated region is slightly larger than the upstream 
direction). The convergence process is very similar for the 
two methods, continuing for about 200 iterations* and even the 
error coefficients are similar. The computation time per itera- 
tion is almost identical for both methods but the memory required 
for MacCormack’ s method is 1.5 times larger than that for the 
Brailovskaya method. 


*It should be mentioned that in the MacCormack scheme (in the 
initial conditions tested) there is one iteration in the x-direc- 
tion for each 8 iterations in the y-direction and the basic At 
in the y-direction is about 10# larger than the At of the 
Brailovskaya method, so that the number of iterations is only 
approximately the same. 
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The fact that the errors in the two schemes are of the 
same order of magnitude (about 0(10- 4 )) while those in 
MacCormack'-s [67] (in which the other components of the 
computation, like boundary conditions, initial conditions and 
also the size of the computation field, are different but 
where the mesh size of the grid is similar) method were of 
the order of magnitude ~0(10- 6 ), strengthens the assumption 
that the main cause for computation errors is inherent in the 
boundary conditions of the external flow (this is also the 
main difference between MacCormack 's [67] original work and the 
MacCormack scheme applied in this study). We will see that the 
distortion in the fulfillment of the flow equations, in the 
region of contact between transition conditions of the incident 
shock wave and conditions of the characteristics, is relatively 
large and the truncation error that it creates (which extends 
inward into the flow field) is the main cause for this difference. 
But, as already mentioned, accuracy of the actual results does 
not fall short of that in the original study by MacCormack [67]. 

It thus remains to choose between the use of the Brailovskaya 
method or the MacCormack method for this research. Prom the 
point of view of accurate results it was preferable to choose 
MacCormack' s method because of its slight advantage, but because 
the computer program was designated for a large scale exercise 
(of about 100 runs) to carry out parametric investigation of the 
results, it was important to use the program with the smallest 
memory that would permit quick and efficient operation from the 
technical point of view. For that reason the Brailovskaya method 
was chosen for the continuous work in this study. 

4.1.2 Comparison of Various Improvement Options 

In this paragraph various options and changes in the basic 
computation method by Brailovskaya, which was chosen for use in 
this study, will be compared. 

The comparison will be made with reference to the experi- 
mental results by Hakkinen et al. [41] and with reference to 
the analytical results by MacCormack [67] and Skoglund, Gay [94]. 

The results of those comparisons are shown in Figs. 6-14 
and relate to the pressure division and friction along the plate, 
to velocity profiles at various cross-sections and also to the 
development of stability coefficients. 

A comparison is made for each option or suggestion of change 
separately, with the reference method applied being the basic 
one (which is defined by a uniform grid, second order differences 
over the entire field, continuous computation of the shock 
wave path, without artificial viscosity and with initial condi- 
tions that already include a schematic form of the field of 
interaction and, finally, with a regular roundoff accuracy to 
6-7 decimals). The run employing the basic reference method is 
indicated in Table 1 as run No. 001. 



*1.1. 2.1 Influence of the Computation Grid Form 
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Two main disadvantages of the non-uniform grid are made 
clear In the comparison between a uniform and a non-uniform grid 
in the x-direction (dense in the middle around the shock wave 
and sparse towards the extremities), which is illustrated in 
Pigs. 6a and 12a by the runs 001 and 006. They are: a tendency 
to develop oscillations in the solution close to the region of 
the incident shock wave and also computation errors that are on 
the average much higher and stem, among other things, from a 
truncation error that is greater in a significant part of the 
computation field (as shown in appendix G'). 

In the non-uniform grid there is a strong pressure gradient 
in the x-direction next to the shock wave (which also accompanies 
the oscillations around it), the significance of it being that 
the shock wave is not sufficiently "spread out" by this method. 

The separated region has a length that is similar to that in the 
uniform grid but it is a little further back in the upstream 
direction of the flow, which provides a better exhibition of 
how far back the influence interaction reaches. The velocity 
profiles for both methods are very similar. 

Summing up, it may be said that at this stage the employment 
of a non-uniform grid is not indicated because, in addition to the 
added confusion that its use causes, it does not offer the 
expected results (unless additional perfections were accomplished). 

*1.1.2. 2 Influence of the Order of Accuracy of Finite 
Differences 

Computation of the fourth order differences (only in the x- 
direction) in that part of the field that is underneath the 
point of entry of the shock wave (between the sonic line and 
the wall) is carried out in run 003 and its results are demon- 
strated for comparison in Fig. 6a. The use of this method was 
originally designated to improve accuracy in the local deri- 
vatives for the regions with strong gradients, however, because 
of the distortion in results generated at the line of contact 
between two computation schemes (of the second and fourth order) 
disturbances develop that express themselves, among other things, 
by oscillations of the separation p w (x) underneath the point 
of entry of the shock wave. On the other hand, a slightly 
larger separated region in the upstream direction of the flow 
is obtained than in the computation using the basic method 
(with second order differences in each field). Once again there 
are almost no differences in the velocity profiles. 

We can then sum up in saying that the incomplete "marriage" 
of the two difference systems in the computation field degrades 
the accuracy of the solution; it is therefore inadvisable to use 
the fourth order differences method for this stage (what's more 
since it requires a memory size and computation time that are 1.5 
times greater than those in the basic method). 
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Note: In the first attempt the fourth order differences 759 

were not used for the whole field since it is not advisable ' 

to compute. with them in the vicinity of the shock wave (for 
fear that a local derivative will be dependent on values from 
both sides of the shock wave, which can be of much greater 
influence in the fourth order differences scheme where each point 
is a function of four neighboring points in the x-direction as 
against two points in the second order differences scheme). An 
additional reason is that away from the shock wave such a compu- 
tation is superfluous, because there are no strong gradients in 
the x-direction, which must be suppressed by the use of this 
method. On top of it all, a fourth order differences computation 
for the entire field will cause a significant increase in memory 
size and required computer time (about 5 times). From all that 
has been said above the conclusion is reached that the use of 
this method is worthwhile only for the subsonic part of the boundary 
layer (which is also the more interesting one with regards to 
computation results). 

4. 1.2. 3 Influence of Artificial Viscosity 

Preliminary computation of about 150 iterations of artificial 
viscosity, which was followed by 200 additional iterations in 
the regular computation (without artificial viscosity) is carried 
out in run 005 , with the results compared to the basic method 
in Fig. 6 a. It appears that artificial viscosity improves the 
results somewhat and, in particular, has a beneficial influence 
for better "spreading" of the gradients around the shock wave and 
on the enlargement of the region affected back of the pressure 
division. The size of the separated region and the velocity 
profiles do not change, however. 

The only drawback in this method is that the coefficients of 
artificial viscosity, CCy and CCx, change with the flow conditions 
and in a manner that cannot be computed in advance, so that the 
optimum values have to be determined in each case separately. Fig- 
ure 7 shows the results of computations for three different 
combinations of artificial viscosity coefficients (one of which 
is the optimum one that is used for comparison in Fig. 6 a). Since 
this optimization involves many tryouts for each case, it is 
not practical when a large number of runs has to be carried out 
under different flow conditions* (in which the Mach and Reynolds /60 
numbers, the intensity of the shock wave, etc., change); for which 
reason artificial viscosity has not been introduced in the rest of 
the runs for this study from here on. Just the same, were there a 
way to develop a simple method for forecasting the value of arti- 
ficial viscosity coefficients as function of the flow conditions, 
it would be worthwhile to use this option to increase the accuracy 
of results. 


*In Fig. 7a the results of three cases for different flow conditions, 
with and without artificial viscosity, are compared where the opti- 
mum coefficients are those chosen from the reference case only 
(pe/po=l- 4 ;Re Xs =3xl0 5 ; Mo=2). The negative influence of this 
artificial viscosity illustrates that the optimum coefficients of 
one case are far from being so in another. 
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4.1. 2. 4 Influence of Shock Wave Computation Method 




In contrast to the basic computation method for which the 
shock wave -is in the center of the computation field (and which 
is also known as the "continuous" method) a method was examined 
in which the shock wave is removed from the computation field 
in a way so that its two sides, up to the sonic line, serve as 
boundary conditions as explained in paragraph 3. 3- 1.5 (this 
method is also known as "discontinuous"). That method was 
employed for runs 007 and 008 (where fourth order: differences 
were also introduced underneath the shock wave). The results 
are shown for comparison in Pig. 6b. 

Comparison of the results brings out that the "discontinuous" ' 
method causes a very strong gradient of pressure in the x-direc- 
tion underneath the shock wave, which is accompanied by some 
slight overshooting downstream of the shock wave. In addition, 
oscillations can be seen in the velocity profile below the shock 
wave. In addition, oscillations can be seen in the velocity 
profile below the shock wave. Oscillations in the solution below 
the shock wave increase when fourth order differences are in- 
troduced into the discontinuous method (in accordance with the 
explanation in paragraph 4. 1.2. 2). In spite of all this there 
is a small advantage to the "discontinuous" method, which is 
expressed in some effect of the pressure to the rear that enlarges 
the separated region. 

One of the main reasons for 
oscillation of the solution near 
the shock wave in the "discon- 
tinuous" method is that during 
the iterative computation of 
the boundary conditions down- 
stream of the shock wave there 
is no possibility for any care- 
ful consideration of the re- 
flected streamlines that return 
in the path of the wave (the 
local intensity of the wave in 
each iteration is computed ac- 
cording to the local ratio of pressures from both sides, the 
pressures being average pressures of the vicinity). No accurate 
local intensity can therefore be obtained for the shock wave and 
for the shape of its curvature into the boundary layer (for which 
reason its location is also not accurate). An additional factor, 
contributing to oscillations, is the form of the boundary on 
both sides of the wave, which is built in the form of steps and 
for which the difference computation is not accurate (this 
phenomenon is also familiar from other sources like Roache [86], 
for instance). 
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This final factor Is notably expressed in a non-uniform 
grid where the step structure in the vicinity of the shock wave 
is even denser and the oscillations of the result even greater, 
until It reaches divergence after 100 iterations (as shown in 
run 009). 
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It should be pointed out, finally, that the "discontinuous /6l 

method" in its present form is not as satisfactory as the 
basic "continuous" method when looking at the results but the 
idea of the removal of the shock wave, around which the strongest 
gradients in the field operate, is worth additional research 
since it points a way for conquering the most difficult ob- 
stacle in this computation - the path of the shock wave. 

4. 1.2. 5 Influence of Initial Conditions 

As has already been mentioned in paragraph 3.3.2, it is 
desirable that the initial conditions not be too far distant 
from the solution, to avoid a potential danger of divergence, 
while on the other hand there is no point in being too meti- 
culous about it so that the general nature of the method is not 
lost. According to the initial conditions formulated in para- 
graph 3-3.2, convergence was obtained to an accuracy of 
0(10- 3 ) - 0(10- 4 ) after 200 iterations (as is indicated in Fig. 

10). For initial conditions that are more removed from the 
solution (like a boundary layer that is undeveloped and uniform 
along its entire length, without appropriate "spreading" of 
the shock wave, etc.) a much slower convergence was obtained, 
on the order of 500 - 1000 iterations, in accordance with the 
distance of the initial conditions from the form of the final 
solution. 

4. 1.2. 6 Influence of the Size of Grid Mesh Dimensions 


To verify considerations for the choice of grid mesh dimen- 
sions (which was made in paragraph 3.4.3 and in appendix D') an 
additional computation, with a grid that is approximately two 
times finer, is carried out (where the number of points is 147 x 
48 as against 76 x 25 for the regular grid) by means of the 
basic method. This fine grid is used in run 04l and the results 
are explained in Figs. 8 and 12b. 

The results obtained show a "spread" for the pressure 
gradient that is just a little bit better but, generally, the 
differences in results from the regular grid and from the fine 
grid are of no significance (and they are even smaller than the 
differences obtained by MacCormack [67] in a similar comparison). 
If there are any big differences at all they are indicative of 
a lack of refinement for getting an accurate solution, but the 
order of magnitude of the differences in this case is quite 
significant, considering the range of accuracy expected from 
this program. 

In the fine grid (which requires 4 times the size of memory 
and 4 times the computer time for iteration) convergence was ob- 
tained after about 400 iterations so that the computation with 
it takes about 8 times longer than computation with the regular 
size grid, without getting the reward of significant improvement 
in results. 
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From results of the computations using the regular grid /62 

and the fine grid, a forecast was made for a number of points 
in the grid for an accurate solution of the difference method 
and from it we obtain an evaluation of the order of magnitude 
for the computation error*. Based on that computation, which 
is found in appendix N, the error in the regular grid is on 
the order of 0(10- 3 ) and in the fine grid it is 4 times smaller. 

4. 1.2. 7 Influence of the Computer “Roundoff Errror on 
Accuracy 

As explained in paragraph 3. 5-3*2, an accuracy of results 
is obtained that is 0(10- 3 - 10- 4 ) so that there is not much 
point in increasing the regular roundoff accuracy of the 
computer (which is between 6 -7 decimals) to double its value 
(to 13' decimals) . To demonstrate that conclusively run 031 
was carried out with double the roundoff accuracy. Ffesults of 
this run are almost completely identical (to the fourth or fifth 
decimal) to those of run 001 (with the regular roundoff accuracy). 
Equally, the error coefficients are nearly alike as shown in 
Fig. 14. 

It should be pointed out that computation with double the 
roundoff accuracy requires double the memory size and 1.5 times 
the computer time used for the regular computation. 

Based on the above, all following computations in the 
program will be carried out with regular roundoff accuracy only. 

4. 1.2. 8 Summary of Comparison of Computation Methods 

Comparison of the various types of computation suggested in 
this chapter with the form of the basic computation discloses 
that actually the basic form, without any attempts for "perfec- 
tions" of any kind, is the most efficient and practical to 
operate for a large series of runs under different flow condi- 
tions . 

In spite of the fact that the other methods will give 
slightly better results in certain regions, they offer disadvan- 
tages that could not be overcome at this stage. In addition, 
these methods complicate the computation scheme and, generally, 
increase the computer time and the required memory size. This 
decision is not a final one, of course, and it is quite possible 
that additional basic research will solve part of the problems 
mentioned to make it possible to combine these methods in a new 
and improved computation scheme. 


-*Bee bottom of page 54 [of original.] 
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We will now sum up briefly the basic computation method, /63 

which will be employed in this study from now on.: : 

a) a uniform computation grid in x and y directions, over 

a computation region of 3 x S 0 x 150 6 0 containing 76 x 25 

grid plate. 

b) roundoff error of 6->7 decimals. 

c) difference method of the second order, using the 
Brailovskaya method for each region computed. 

d) no expression of artificial viscosity added to the 
difference equation. 

e) "continuous" computation of the shock wave path, i.e., 
the shock wave in the center of the computation field and 
defined as external compulsion only for the boundary conditions. 

f) initial conditions include an approximate and very 
general description of the interaction field. 

A convergent* solution is obtained with this method (charac- 
terized by asymptotic values for the stability coefficients) 
after about 200 iterations. 

Computer IBM 370/168, on which the runs required for this 
investigation were performed, needs 4 minutes of CPU per run and 
a memory size of 256 K (for continuous computation of the stability 
coefficients about double the time and memory size will be required). 

4,2 Parametrical Analysis of the Results /64 

Analysis of the results of computation for various flow 
conditions, as obtained through the basic computation method, 
will be carried out below. 


*The description of development of the stability coefficients in 
Fig. 10 shows that after about 200 iterations they approach an 
approximate asymptotic value, the average order of magnitude for 
errors of the most frequent variables is about 0(10“3) and the 
order of magnitude for the average derivatives is 0(10“ 2 ). Those 
values are in accord with evaluation of errors as explained in 
appendix L, which is 0 ( 10“3 -lo - ^) of the actual values of the 
variables . 


A description of the division of errors in the field, in 
Fig. 11, acknowledges that the main errors of computation are 
concentrated around the shock wave in the center of the field. 
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*1.2.1 Influence of Plow Conditions on the Interaction 


For most of the flow conditions for which computations 
were made neither experimental nor theoretical data exist for 
comparison; therefore the analysis will be made with the 
assumption that the results are correct, which is based on a 
comparison with the few results available until now. 

The runs described in this part are concentrated in Table 2. 

*1.2. 1.1. Comparison with Previous Experimental and 
Computational Results 


Runs 101, 102 and 103 were made for comparison with exis- 
ting results for the following conditions: 

Mo = 2, Re Xs = 3 x 10 5 , x s = 5 cm 
P e /p o = 1.2, 1.-4, 1.9 

^Comparison is made with experimental results of Hakkinen 
et al. [ 4 1 ] * for all three cases and also for computation results 
of MacCormack[67]** and Skoglund, Gay [9**] for p /p =1.2, 
l.*i only. e ° 

In Pig. 15 a comparison is made for p e /p Q = 1.2 where the 
flow has not yet separated. The results of run 101 are very 
close to those of .'MacCormack[67] (except for the slightly 
stronger pressure gradient and the slightly greater friction 
downstream, in the computational results). 

The most important comparison was made in Fig. 17 for p /p 
= 1.**, where the flow has already separated. It turns out tHat° 
for all computational methods (including the present one 
demonstrated by run 102) a lower pressure is obtained upstream, 
when compared with experimental data***, just as the boundary 


*From a description of the experiment and attached Schlieren photo- 
graphs it is certain that the experiment was performed under 
laminary flow conditions. 

**Results of MacCormack’s computations in this paragraph are from 
his basic study, MacCormack[ 67 ] , as contrasted with results of 
the present investigation in *1.1.1, using MacCormack's scheme. 

***The difference arises perhaps from the inaccuracy of measurement 
in the separated region (another possibility is a computational 
error common to all methods of computation but that does not 
appear likely in view of the successful solutions for a variety 
of problems, obtained through these computation schemes). 



layer is smaller below the shock wave.* 


The difference between results from run 102 and the 
other computational data is not a big one. In consideration of 
these small differences it may be said that the MacCormack 
method is the best one, that of Skoglund, Gay [94] is the least 
good one and the computational method used here is somewhere in 
the middle. 

In general we will see that the present method shows a 
slightly larger pressure gradient below the shock wave and 
less pressure increase in the upstream flow, as compared to 
MacCormack* s[ 67 ] method, but a more accurate pressure (compared 
to experimental data) for the downstream flow after the shock 
wave. In the present method the separated region is slightly 
smaller and the friction coefficient for the downstream flow 
is greater. There are no significant differences in the velocity 
profiles of all three methods of computation. 

Fig. 19 explains the comparison of data from run 103 for 
Pe/Po = 1«9 to experimental results alone (computational results 
are lacking). The pressure division is perfectly correct 
but even here the upstream pressure increase is smaller and the 
separated region smaller, when compared with experimental data. 

Friction of the downstream flow is also lower. As with results 
for p e /p Q =1.4 the velocity profile below the shock wave, as 
obtained through computation, shows a thinner boundary layer 
and a separated region that extends less into the field. 

In addition to the comparisons listed above, flow field 
charts with iso-Mach lines, iso-flow lines and iso-pressure 
lines were plotted in Figs. 16, 18 and 20, for runs 101, 102 
and 103. These charts present a good picture of the flow field, SLL- 

which compliments the obtained results from Figs. 15, 17 and 19. 766 

The charts clearly show the separated region, the incident shock 
wave and the reflected streamlines, as well as the development 
of the pressure change in the field. 


* This phenomenon must also be laid to an erroneous measurement 
in that region. There is an additional possibility, namely 
that during the experiment the flow in this region was subject to 
characteristics of turbulence that changed the y gradient of 
the velocity profile, leading to an increase of the separated 
region (while all numerical computations are only for laminar 
flow). As a matter of fact, a transfer of the experimental 
velocity profile from Fig. 17, in the y-direction, up to the 
coalescence of the points where u = 0, shows a much better fit 
between the experimental profile and those obtained through the 
computation methods. 
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To summarize, it has been shown that a comparison of data 
from the present method of computation, with existing results 
from computations and experiments, shows them to be entirely 
reasonable and falling only slightly short of those from the 
competing MacCormack [67] method (in the range of data where 
such a comparsion was possible). For strong interaction, for 
which only experimental data were available for comparison, it 
appears that the separated region is smaller but the form of its 
development, with the intensity of the shock wave and the Mach 
number, is in general accordance with experimental data (as we 
shall see during the continuation of this data analysis). 

The advantage of the present method is in its relatively 
early convergence (apparently with a large error of computation, 
which does not, however, degrade the accuracy of the final result 
very much), which permits its practical application for large 
numbers of runs with different flow conditions. 

Note: as already mentioned in paragraph 4. 1.2.3* some 

slight improvement in the accuracy of results might have been 
achieved by the use of artificial viscosity, but, for the 
sake of simplicity of operation and aviodance of additional 
computations to find its optimum coefficients, it was not done. 

In the following paragraphs the parametric analysis of the 
influence of the various flow conditions on the interaction 
field will be performed (through the runs illustrated in Table 2). 

4. 2. 1.2. Influence of Shock Wave Intensity 

Fig. 21 describes the influence of the shock wave (in the 
range of pressure ratios 1.2 <_ pe/po £ 3.15 s ) on the inter- 
action when Mo = 2; Rex s = 3xl0 5 (in runs 101, 102, 103 and 104). 

The increase in intensity of the shock wave enlarges the region 
to the rear that is influenced by the pressure and, along with 
it, the size of the separated region also increases towards the 
upstream direction to the flow. At the same time the separated 
region enters deeper into the field and because of it the edge /67 

of the boundary layer moves further away from the wall. Fig. 

22 describes the flow field picture in run 104 where a very intense 
shock wave of p e /p Q = 3-15 exists. We can see the change in 
flow field characteristics more clearly by comparing it to 
weaker interactions (as in Figs. 16, 18 and 20.) 


*In this comparison no shock intensities higher than for Mo=2 
were included and for p e /Po = 3-7 a subsonic region starts after 
the shock wave (run 105) and the boundary conditions of the 
characteristics cannot be applied. 
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It must be stated that in reality a transition of the 
boundary layer may develop so that it will be turbulent for 
intense shock waves (particularly downstream of the shock wave) 
so that the (laminar) computation results won't be so accurate; 
but the main interaction phenomena take place upstream of the 
flow and are not much influenced by that . 

4. 2. 1.3. Mach Number Influence 


Pig. 23 describes the influence of the Mach number on 
interaction in the range* 2 <_Mo <_ 4.5, when p e /Po = 3*15 and 
Re Xg = 3 x 10 s (runs 104, 111, 112 and 113). The comparison 
was made for a relatively intense shock so as to obtain the 
Mach number influence at strong interactions. 

It appears that the higher the Mach number the more acute 
the angle at which the shock wave enters the field, which weakens 
the interaction. That is visible from the results, inasmuch as 
the higher Mo the more the region to the rear, subject to influence 
of pressure and of size of the separated region, will decrease. 

Fig. 24 describes the flow field of run 112, with a high 
Mach number of Mo = 4. The weakening of the force of interac- 
tion, when compared to lower Mach numbers (as in Fig. 22 for 
instance), is clearly visible in this picture. 

4. 2. 1.4 Reynolds Number Influence 


Fig. 25 describes the Reynolds number influence on interac- 
tion in the range 1 x 10 4 <_ Re Xs £ 1 x 10 6 , when Mo = 2 , 
p e /p 0 =1.4 (runs 121, 122, 102 and 123). 

The results show that a decrease in the Reynolds number 
pushes the pressure increase downstream because of the shock 
wave. This is in agreement with a similar investigation carried 
out by Reyhner, Fldgge-Lotz [84] for Mo = 3. The separated region 
increased with the change in Re x from 1 x 10 4 to 3 x 10 5 but de- 
creased gradually when the Reynolds number was further increased. 
One possible explanation for that is that the computation method 
tends to show (though in a very coarse manner) the solution for 
turbulent flow for this range of the Reynolds number, while in 
reality a transition towards turbulent flow in the boundary 
layer exists for this range (around Re Xs “lx 10 6 ), in which 
changes are much smaller because of the much stronger flow 
resistance to pressure gradients (which also decreases the in- 
fluence of interaction and, with it, the size of the separated 
region) . 


*For higher Mach numbers (for the same pressure ratio p e /p Q = 3.15) 
instabilities start to develop in the solution. This is expressed 
in run 114, for p e /Po = 3-15, Mo = 5, where the solution diverges 
after 150 iterations. 


61 



Together with the above, the difficulties of exactly de- 
termining the size of the separated region from computed results 
must be pointed out. This is due to a somewhat unclear de- 
termination of the border between the separated region and 
the boundary layer (because the changes in velocity components 
are very small for this region). 

2, Pig. 26 describes the flow field in run 121 for Re x = 1 x ° 
10 , showing the influence of a low Reynolds number on the inter- 
action field as compared to the influence of higher values (as 
in Fig. 16, for instance). 

4. 2. 1.5 Influence of Shock Wave Entrance Location 


Fig. 27 describes the influence of a change in the location 
of shock wave impingement for the range 3 cm < x s < 7 cm, when 
Pe/Po =1.4 and Mo = 2 and the Reynolds number of the reference 
point x. = 5 cm is 3 x 10 5 . The results were obtained from 
runs 131, 102 and 132. 

It becomes clear that the larger the distance between the 
impingement location of the shock wave on the plate and the leading 
edge, the smaller the pressure gradient along the plate and the 
greater the influence of pressure towards the rear and on the 
size of the separated region. In fact, this is the result of 
local growth of the boundary layer thickness in the interaction 
region with the increase in x s , so that the boundary layer 
becomes more sensitive to pressure gradients. 

4. 2. 1.6 Prandtl Number Influence 


Fig. 28 describes the comparison between results for Pr=0.72 (val- 
ue for air) and 1.0 (an approximate value used often because it 
simplifies part of the equations), when Mo = 2 , p e /p 0 =1.4 and 
Re Xg = 3 x 10 5 '(runs 102, 141). It turns out that neither 
pressure division, the friction coefficient, or the separated 
region show any change because of it. The only difference ex- /6Q 

pressed is a slight change in the temperature profile next to 
the plate, which causes an adiabatic wall temperature that is 
about 5 % higher for Pr = 1.0. It is to be expected that this 
influence will grow with the Mach number and that, because of 
it, the form of the boundary layer and the friction coefficient 
will change slightly. 

Those results were found to be in good agreement with a 
similar study by Reynher, Fltlgge-Lotz [84] for Mo = 3. 

From this investigation it may be concluded that the ap- 
proximations of Pr = 1.0, that are occasionally made in various 
studies, are generally justified and particularly when the 
Mach numbers are low. 
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4. 2. 1.7 Summary of Flow Condition Effects 




From examination of the parameters in this section it 
appears that the pressure ratio (i.e. the shock wave intensity 
exerts the most important influence on the interaction, es- 
pecially on the size of the separated region. 

The Reynolds number, too, influences results considerably, 
particularly since the width of the boundary layer and its 
sensitivity to pressure gradients depend on it. 

The other parameters like Mach number, location of the 
shock wave impingement and the Prandtl number, exert much less 
influence if at all. 


Fig. 29 describes the influence of parameters Mo, p e /p Q and 
Re Xs on the development of the boundary layer along the interac- 
tion field. As expected, the Reynolds number is the primary 
determinant for shape and width of the boundary layer develop- 
ment while the influence of p. /p and Mo is of a much smaller 
order. e ° 


4.2.2 Influence of Boundary Conditions next to the Plate 770 

(Heat and Mass Transfer) 

r The influences of heat transfer (maintenance of the plate 

at a fixed non-adiabatic temperature) and mass transfer (suction 
or self-injection into the plate) on the flow field in general 
and on the separated region in particular, when all other flow 
conditions (only in this part) remain fixed at the values Mo=2; 
p e /p 0 = 1.4; Rex s - 3 x 10 5 , will be described further on. The 
runs for this section are described in Table 3. 


4. 2. 2.1. Influence of Heat Transfer Next to the Plate 


The influence of maintaining the plate at various fixed 
temperatures in the range* 0.4< T, f < 1.4 (when T^=CpT w #/U 0 # 2 ) , 
as obtained from runs 201 and 206j is described in Fig. 30. 


Maintaining the plate at a fixed non-adiabatic temperature 
causes a heat transfer next to the plate that influences the 
temperature profile but, as it appears from the drawing, the 
influence is localized only and attenuates quickly until it 
disappears entirely at about the middle of the boundary layer. 
The thermal conductivity coefficient next to the plate changes 
according to the plate temperature; its value is negative for 
a plate temperature that is higher than the adiabatic one and, 
conversely, positive when it is lower. 


*The temperature of the adiabatic wall for Mo = 2 is T =1.05. 

w 
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On the other hand, the pressure division, the boundary 
layer development and the size of the separated region as result 
of temperature changes along the plate, are hardly influenced 
at all. 

Those results are in agreement with data computed by 
Reyhner, Fltigge-Lotz[84] for interaction that ends next to a 
blunt corner. 

It may be concluded that there is no practical way to 
control the size of the separated region through heat transfer 
along the plate. 

4.2.2 Influence of Mass Transfer Close to the Plate 

In contrast to heat transfer, mass transfer close to the 
plate (which is expressed here by introduction of suction, or 
injection, normal to the plate length) does have significant 
influence on the flow field of interaction, and particularly on 
the size of the separated region*. 

Fig. 31 describes the influence of velocity, normal to 
the plate and next to it, in the range of - 0.03 < V * + 0.02 

41 w 

(when V = V ) , as obtained from runs 211 to 2i5. 



What stands out in particular is the strong influence on 
the friction coefficient, the size of the separated region and 
the profiles of the boundary layer. The suction next to the 
plate causes the boundary layer to "stick" to it by reducing the 
separated region until it disappears completely (in this case 
at V w = -0.02). Injection, on the other hand, causes a broadening 
of the boundary layer and an increase in the separated region. 

The friction coefficient increased with the increase in suction 
velocity and decreased with the increase of injection velocity 
into the plate. 

Pressure division is almost unaffected, except for slight 
oscillations in the initial and final regions of that section 
of the plate where suction or injection are present (these 
oscillations, which are not visible in Fig. 31, are apparently 
the result of local numerical discontinuities in the boundary 
conditions on the plate). In addition to what was said above, 
suction influences the slight delay in pressure increase down- 
stream of the shock wave while injection causes slight accel- 
eration upstream of the shock wave. 


^Parametrical analysis of the influence of suction on the 
various flow conditions and the conditions required to prevent 
separation of flow are shown in the following section 4.2.3. 



Fig. 32 describes the flow field influenced by a suction of 
-0.02 (-run 212). The figure shows clearly how the flow >J 

lines bend and the separated region disappears, when compared 
to flow without suction influence under the same conditions (as 
in Fig. 16). 

Fig. 33 describes the influence of suction and injection on 
development of the boundary layer shape. 

Finally, about influence of boundary conditions on the 
interaction it can be stated that heat transfer has almost none, 
while mass transfer has very significant influence on the size 
of the separated region and on the shape of the boundary layer. 

4.2 .3 Detailed Analysis of Suction Effect on the Separated /72 
Flow Region 

This section goes into more detail to explain the influence 
of suction on the separated region for various flow conditions 
and the amount of suction required to prevent separation for 
those conditions. The subject receives a great deal of attention 
in comparison with other parametrical studies, since it is of 
important practical significance. 

4.2. 3.1. Influence of the Location of Suction Along the 
Plate 


In runs 211 to 229 the influence of changing the location 
where suction on the plate starts in the region* 2.6 cm <_ x <_ 4 . 4 
cm (when termination of suction is always -0.03 >V W > 0, with 
flow conditions being: Mo = 2 , p e /p Q = 1.4 and Re Xs = 3 x 10 5 . 

The results show that this has no actual influence on the solution 
and they are nearly identical for all cases. 

Since in this comparison the suction region always includes 
the entire separated region, the assumption was made that only 
suction in the separated region proper had any influence on the 
results (mainly because of the size of the separated region), 
while no change at all was to be expected in the results due to 
suction from any section outside it. 

To verify this assumption, an additional series of computer 
runs was arranged (runs 251 to 267 described in Table 4) for 
Mo = 2, p e /p 0 = 1.4, Re X £ = 3 x 10 s and V w = -0.02. For those 
runs the start and finish locations of suction changed over a 


*The suction region always starts before the start of flow separa- 
tion along the plate, in this case. 
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very broad range, so that many cases could be Included where 
only part of the separated region was Influenced by suction 
while a part was not. The data shown in Figs. 34 and 34a 
confirm that suction is effective (as far as its influence on 
the separated region is concerned) only when it operates 
within the region but not outside it. To benefit most from 
suction it must therefore be maintained over the entire 
separated region. 

Since it is not known in every case just where the separated 
region will appear, the following parametrical tests were run 
for suction. The region of suction is thus very broad and covers 
most of the plate length so as to obtain maximum benefit from 
its influence for each case. 

4.2. 3.2 Influence of Suction for Various Flow Conditions 


In the series of runs described in Table 5 (runs 301 to 
452) the influence of the following parameters was examined: 

5 x 10 4 < Re Xs £ 1 x 10 6 ; 2 <_ Mo < 4; 1.2 £ p /p Q 1 3.2 and 
3 cm x s >_ 7 cm^ with suction' changing over the range - 0.3 > 
V w >0. Results of the runs are shown in Figs. 35, 36 and 37, 
with the size of the separated region as a function of suction 
velocity for various flow conditions and in non-dimensional 
form. 


The figures disclose that the pressure ratio Pe/Po influ- 
ences the increase of the separated region and so contributes 
to the increase in suction velocity required to eliminate it. 

It should be noted here once more, as has already been said in 
paragraph 4. 2 . 1.1 that the separated region obtained through 
computation is smaller than the one in reality, particular for 
high pressure ratios; however, the general purpose of the in- 
fluence exerted by the pressure ratios on the separated region, 
for various suction velocities, appears to be superficial and 
reasonable . 

The influence of the Mach number tends toward decrease of 
the separated region, so the suction velocity required to elim- 
inate it decreases as the Mach number increases. 

The influence of the Reynolds number is more complex. As 
already described in paragraph 4.2.1. 4, the separated region 
starts to decrease when Re Xs exceeds the value of 3 x 1CP 
(while in the lower range the separated region increased together 
with Rex s ). On the other hand, when the suction intensity 
increases (V w < . 18 ) the opposite situation is created, with 
the separated region continuing to increase with the Reynolds 
number, at least to Rex s = 10 6 . If we base ourselves on the 
estimate offered in paragraph 4. 2 . 1.4 we can claim that an 
increase in suction causes rejection of the tendency for 
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transition and with it for turbulent flow, so that suction is 
responsible for the return to laminar flow, which is apparently 
also expressed in the numerical solution. To sum it up, the 
more the Reynolds number goes up the stronger is the suction 
required to eliminate the separated region. 

The impingement point of the shock wave exerts next to no 
influence. In the figures the non-dimensional results for x s 
vary as they combine, with good approximation, to one line. 

(That result is forced by circumstances since the length of 
the non-dimensional separated region is also referenced to x s ). 

4. 2. 3. 3 Suction Required to Eliminate the Separated 774 

Region 

Figs. 38, 39 and *10 (which were constructed from the 
previous Figs. 35, 36 and 37 ), describe the dependence of 
suction velocity, needed for elimination of the separated 
region, on each of the main parameters (pressure ratio, Mach 
number and Reynolds number). 

With the aid of appropriate extrapolation from those 
figures a summary table was drawn up in Fig. *11 where the 
needed suction velocity for prevention of separation is given 
as the function of those three parameters. 

It is also appropriate to recall that Carter [13], in his 
investigation of the influence of suction on the flow around a 
rounded corner at an angle of 10 ° (for which he got a pressure 
ratio of p e /p Q = 2.2) with Mo =3 and Re Xg = 1.68 x 10 4 , found 
a requirement for suction estimated at = 0.1 to prevent 
separation; despite the difference in the form of interaction 
in this research, it turns out that this value is very good 
agreement with that received from the table in Fig. *11 for such 
flow conditions. 

*1.2.3.*! Approximate Formulation of the Suction Needed 
To Prevent Separation 

In the computation described in appendix 0, a single- 
valued connection was found between the characteristics of the 
separated region and the suction velocity needed for its 
elimination. 

The separated region is characterized by its dimensions 

Ax s ; Ay s and by the maximum negative value of the flow function 

(with the focus of the bubble in its center) wich is 7 ' 

mx. 


The connection is, to a good approximation. 





- t+. 


(4-1) 
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where the fixed value is in the range 0.3 t 0.^2, the bubble 
area is defined as v = Ax g . Ay s .1 and the flow function 
7 = v/p Formula (4-1) was tested for 10 different combinations 
of flow conditions (Mo;Re. v ; p /p ). 

x s e ° 

5. Discussion and Conclusions /75 

The finite difference computation method, investigated 
in this study for solving problems of interaction between . 
shock wave and laminar boundary layer (through solution of the 
complete Navier-Stokes equations), provided excellent solutions 
when compared with the experimental and theoretical information 
presently available. The simplicity of the method's operation 
and the relatively short time needed for the solution, permitted 
the execution of a large number of runs for various flow condi- 
tions, from which it was possible to learn about the interaction 
characteristics and the principal factors that influence it. 

Special emphasis was placed on researching the influence of 
suction on the prevention of flow separation next to the plate, 
which is of great practical significance. 

We will now discuss results and main conclusions reached 
from this study: 

f" 5 . 1 Computation Method 

The method of computation demonstrated in this study is 
characterized by the following factors: 

a) The method is constructed simultaneously with the second 
order difference schemes of Brailovskaya and MacCormack, for a 
uniform grid. 

b) The field of computation stretches over a range of 1506o 

length and 36 0 width around a region where the shock wave impinges 

on the plate. The field is divided into a grid of 76 x 25 points 

that forms uniform meshes with the dimensions 26o x 0 . 125 . 

Convergence for a normal computation occurs within about 200 
iterations and the obtained accuracy for the variables (through 

truncation errdr) is in the range of 0(10- 2 3 ) -0(10- 4 ). Such a 

computation, carried out on an IBM 370/168 computer, took 4 
minutes of CPU time. 

c) Boundary conditions for this computation are: 

1. Inlet cross-section - profile of boundary layer 
per Polhausen. 

2. Exit cross-section - zero gradients in the flow 

direction. 
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3. Plate - no velocity components (non-slip condition) 
and no gradients normal to the direction of temper- 
ature (adiabatic plate) and pressure (for approxi- 
mate computation of density). 

4. External flow - according to the characteristics 
from the field (except for the two sides of the 
shock wave), this boundary condition distinguishes 
the present computation method, with reference 

to the other existing solutions, and it permits 
the reduction of the computation field up to near 
the boundary layer. 

d) The computation results agree in the main with those 776 

known from experiments, except that a much smaller separated 
region is obtained for very strong interactions, but the tendency 
in the development of this section with changing flow conditions 
is in the right direction and the right proportions. 

The results are in good agreement with those from previous 
computations and fall only slightly short of MacCormack's results, 
while showing improvement over those by Skoglund’, Gay (those 
two studies were solved with schemes similar to the one used 
here, though the external boundary was assumed to be far from 
the plate and fixed conditions were given for its entire length, 
from both sides of the shock wave. The solution in those studies 
converged after several thousand iterations to a truncation 
error of less than ~0(10~ 3 * * 6 ). 


e) A number of additional auxiliary methods were tried 
during this study, but were not introduced into the final method 
of computation because they have not yet shown themselves suitable 
for that. They are, however, worth additional research so that 
they may become helpful in improving the results. 

1. Artificial viscosity - an option was prepared for 
addition of artificial viscosity terms to the 
difference scheme, which improves the results 
slightly but requires computation for optimiza- 
tion of the numerical coefficients. 

2. Non-uniform grid - with changing dimensions in 
the x-direction so that it is dense in the center 
(next to the shock wave) and sparse towards the 
ends. The solution for such a grid has a tendency 
for oscillation near the shock wave. 

3. Fourth order finite differences - these were used 

for the x-direction, below the shock wave, but 

caused oscillations in the zone of contact with the 

region of second order computation (in the vicinity 
of the sonic line). 
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Removal of the shock wave from computation field - 
this "discontinuous" computation of the shock 
wave by fixing its two sides as boundaries of the 
field, did not bring the expected results and 
caused oscillations around the shock wave. 

5.2 Parametrical Analysis of the Interaction /77 

The influence of shock wave intensity, 1.2 <_ p /p <_ 3.2 and 

its location of operation 3 cm < x 7 cm, were investigated. The 

s < 1 x 10 6 

Mach numbers 2 <_ Mo <_ 4.5, Reynolds numbers 1 X10 <_ Re 

A 

S 

the Prandtl numbers Pr = 0.72, 1.0 were used. Also investi- 
gated was the influence of heat transfer in the range 0.4 < T 

_ 1.4 and mass transfer in the range -0.03 _ V _ 0.02 along 
the length of the plate (non-dimensional values^. 

The principal results obtained were: 

a) Increase of pressure ratio in the flow field enlarges 
the separated region and the pressure influence to the rear 
(pressure increase in front of incident shock wave location). 

b) The further the location of the incident shock wave 
from the leading edge the more "spreads" the influence of the 
shock, which is expressed in the decrease of the pressure gradient 
and enlargement of the separated region. 

c) Increase of the Mach number weakens the force of the 
interaction and reduces the separated region. 

d) Increase of the Reynolds number for the point of inci- 
dence of the shock wave enlarges the separated region, but above 
the value of 3 x 10 5 it reduces it. Also, the pressure gradient 
is moved forward in the upstream direction of the flow. (Because 
of the possibility that at high Reynolds numbers the flow may 

no longer be laminar past the shock wave, it is not certain 
whether the results for this region are reliable). 

e) Influence of the Prandtl number on the result is 

nearly negligible, except for a small local change in the tempera- 
ture profile in the lower part of the boundary layer. 

f) Mass transfer (suction or injection) next to the plate 
influences the size of the separated region greatly; strong 
suction can completely elminate the separation and reduce the 
boundary layer thickness significantly while injection brings 
about the opposite results. Pressure division and the rest of the 
variables in the flow field are almost unaffected. 


70 



m. 




Because there is great practical implication to the reduc- 
tion of the separated region during the interaction, many suction 
conditions were investigated in various combinations with the 
flow conditions so as to gain detailed information about. its 
influence. 

First it was found that the suction influence is effective 
only in the separated region proper; active suction along this 
entire region must therefore be aspired to. 

It turned out that the suction values required for preven- 
tion of separation are in the range 0> V > - 0.04, for the 

w 

range of various flow conditions investigated (per the previous 

paragraph). A detailed table and charts were prepared of the 

suction velocity required for prevention of separation as 

function of the main flow characteristics Mo, Re and p /p ; 

’ x e o’ 

s 

also an approximate empirical formula was found that connects 
the above suction velocity to the characteristics of the separa- 
ted region (its dimensions and the value of the flow function 
within it). 

It is to be noted that, in spite of inaccuracies in the 
computation of the separated region during the strong interac- 
tions, the results are methodical and consistent and so permit 
examination of the relative influence of the various suction 
conditions on the size of the separated region, for various flow 
conditions . 

At any rate, it is clear that one cannot rely on the results 
(the suction values needed to prevent separation) with complete 
confidence, but only as a relative size and as approximate order 
of magnitude. 

The two principal results in this paragraph - effectivity 
of suction within the confines of the separated region alone and 
the order of magnitude of the suction velocity needed 'to prevent 
separation, which is up to 4$ of the expected speed - are of 
important practical significance for the design of aircraft 
components where separation of the boundary layer must be prevented 
during interactions of this kind (particularly at jet engine 
inlets and at control surfaces). 



From the results of this research, as well as from similar 
studies carried out in the past or recently, it appears that it 
will be worthwhile in the future to concentrate on the following 
subj ect s : 
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a) Improvement of those numerical methods that have not 
yet been developed sufficiently (as mentioned in paragraph 5.1), 
so as to increase the accuracy of their results. 

b) Extension of the computation method to turbulent flow by 
introducing the appropriate model of turbulence. This subject 

is particularly important because turbulent flow occurs in a 
significant part of practical interaction problems. 

c) Formulation of the computation method into a general 
form so that it would be suitable to a wide spectrum of interac- 
tion problems like a compression corner, forward and rear steps, 
etc. The method's utility should be such as to lend itself for 
quick and easy use for all cases. 

d) Adaptation of the computation method to any combina- 
tion of geometric boundary conditions so as to permit investiga- 
tion of flow over complex bodies on which a number of forms of 
interactions of different types occur (so the method will become 
a sort of "numerical wind tunnel"). 


This final subject will not be realized in the near future 
since it requires a lot of development and the use of a memory 
size and of computer time that are so large as to be impractical 
at this stage. But with future development in computers, it is 
to be expected that the matter will then be more practical. 


Appendix A* - Non-dimensional definitions of characteristic flow 
parameters 


a. Mach Number 
By definition 


M - 




( A-l ) 


by setting 

- U.VA.0 

V* “ V lA-o 

T* --- T U-o V C P 

= Cp-Cv 
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we get 


( A-2 ) 



b. Plow function 


We will define a non-dimensional flow function through 
non-dimentional variables below: 

f -- - f/ 

I J 


Mr ^ _f 


(A-3) 


and we will get the change of the flow function between two 
points through the linear integral 

(y ft 

Wo - rw t ' /ix) -f-j. J (a-h) 

When is the average density between points a and b 

c. Friction Coefficient 
By definition 


CP - 


«/4 JV v.t 


'C * > K 


(A-5) 


by setting 






we get 




U.* = \X.VA. o! 
= JA jAo 


(A-6) 


d. Heat transfer coefficient 


By definition 




f f v.-T To 


(A-7) 


* k* 

M * 
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By setting 


It i 

Po V.tT M S 


“ Pr 


we 


Eet wb:. ^'4(i?y 


= y> jKo 

t*~ r 

tf*'}** 


with the reference condition V = 0 and so joj 


+ v . 


the following is valid 
from which we get 


,«.= Cf- <v X- pt/c* 
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Pr 
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A ■*' 


u* 


2ff' T c 


(A-8) 


/8l 


7 ^ 



Appendix B 1 - Flow profile computation of inlet section (per 
Polhausen method) 


Assumptions for the computation are: 


a) 

*x± = ° 


(B-l ) 

b) 

fir*) 


(B-2 ) 

c) The relation between the Drofiles of temnerature and 
velocjtv in the laminar boundary laver and comoression on a 
f]at Plate 3 s (ner Schlichtlng fQp.l'): - 


T . , r— y~! 

r^ r To ~ .z 

* [<- (£-)"] - £ 

j (B-3) 

We 

define the momentum equation 

as 



cLo* /jl w f y ^ \ (b-4) 


when momentum thickness is given by 


also 


£■*(»•-) . „ , 
»(»•>* / -fe- tzj* r *.*) * 

“ rU*’) 

£V, 


(B-5) 


(B-6) 


We define a new variable 


i- v «? I 


J" 


and get * 

'S 


--if 


d« 


•’■Wf 


Through exchange of variables from y* to n 
we get (B-5) again 


/ M* e <> ^ 


(B-7) 


dr, ( B— 8 ) 


(B-9) 







and also valid is 


\i£/ U* = K.'Vfe,, « *.*£ /*<*», 


(B-19) 


From substitution of (B-3) and (B-12) into (B-8a) and 
integration over y# we finally get the following non-dimen- 
sional expression 




(B-20 ) 


Formulas (B-3), (B-12) and (B-20) define the flow profile. 

As to component V of the velocity at the inlet section V = 0 
can be assumed, based on the following evaluation: 

V (y) fulfills the order of magnitude ’ # . nP (B-21) 

° * V I ) * U. — 

We then calculate the numerical value of this order of 
magnitude through the use of the following data (which 
demonstrate the computation condition in this study) 


he - o>. R<c % * 

we get 

>i s *» 


V- 


_ s ’ (V-i)Hc' T /... A ,, i._ r*l‘- ^ ■ 

jrvy- L- h - • *j » 

sal 

through substitution in (B-20) we get 

{(&){■&)(■&>) [•’- sr'"* 

j A as O.olSi- 


tfe -o.osrm. .y^ -o.oiT 


( B-22 ) 


(B-23) 




V (B-2H) 


and from (B-15) we get 




T>X i*, 

from the condition at the end of the boundary layer 
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is valid 


\f_ i£ 

\ j . = <j a where AX~* 
we get 


c < 




(B-25) 


(Note: in the region of flow conditions 2< Mo < 5 and 
10 4 <Re < 10 6 there is no change in the order of magnitude of 
V(y)). 


Prom this evaluation it becomes clear that V may he neglec- 
ted at the inlet section with good approximation and we may assume 
that 

V = 0 (B-26) 
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Appendix C* - 
pressure ratio 


Evaluation of shock wave force for the required 
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The pressure ratio parameter P e /P Q characteristic for 

interaction (next to parameters x g . Re , Mo) expresses the 

s 

pressure ratio between exit cross-section and inlet cross- 
section. This pressure ratio is generated by. the transition 
along the incident shock wave and the diffuse reflected flow 
lines (which combine into the reflected shock at great distance 
from the plate). As a practical matter, either experimentally 
or theoretically, the pressure- ratio P 2 /P-j_ from both sides of 

the incident shock wave at the external boundary of the field 
is furnished and as result , the ratio P e /P 0 is obtained. 

It is generally not possible to calculate P 2 /Pq directly 
from P e /P Q and so an iterative computing procedure is required. 


We will show a sample procedure for a simple case dealing 
with a shock wave that touches the plate and is reflected from 
it (without consideration of the boundary layer). 


G 


iven data: p^/p^, M 1 


incident wave 


reflected wave 


Required: P2 //p i 


Solution procedure: 



a. an initial value is assumed for p_/p , 
for instance 


b. from the formulas for an 
oblique shock wave we get 
the angle of the incident 
wave 


-i • -* 
0 t ~ 


.nr 


Hi 


c. angle of 
after the 


flow deflection £ 9i 

shock 


J 


d. angle of the reflected wave 


Oz* Q, 


-nl 


e\ Mach number downstream 
of shock wave 


am ^ r( H.^9.r-+7rT 

*■' S-16,4)! tX 4] 


f. Pinal pressure ratio 


(C-l) 
( C-2 ) 
(C-3) 

(C-4) 

(C-5) 

( c— 6 ) 


when 


i^v - — • - - — ^ 4}] 

Veil “ <\ f!. p ‘ L 
Np) i 


> £. 


then (P 2 /Pq) is corrected through 
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1 

e* s i 

I 4 

i p k 

A, 

(%•> 

T. + £ )' 

% 

p 1 

and from here we return to stage b. ' 1 

1 -- 

a 


£-0 


This process converges quickly and it becomes clear that 
the first guess in a. is very close to the final solution. 


When speaking of the region close to the plate (not far 
from the edge of the boundary layer) there is as yet no re- 
flected wave but only a series of diffuse compression expansion 
flow lines, as explained in paragraph 1.2. 

Since the form of the flow lines and their direction are 
not known in advance, it is impossible to set up a computation 
like that described for the model. 

Therefore, it is advisable to assume an initial pressure 
ratio- p^/p 1 (an approximation of p 2 /p 1 = * / f^7Pq" is recommended) 

and to find the resulting (p./p.), from the flow field solution 
and then to correct P 2 /p 1 based'on the amount of divergence of 

(p^/p^), from the given value of (P /p^). 

It must be pointed out that it is more important (particu- 
larly for the purpose of comparison with other results) to 
obtain an accurate p /p ratio (indicated in the sample by 

/ \ 0 O 

p 3 / ^l^ than the right combination of waves and diffuse flow lines 
that generate this value. 

An additional factor 
influencing p 2 /p 1 is 

the nature of the inter- 
section between the in- 
cident shock wave and 
the reflected diffuse 
flow lines, determined 
by the distance between 
the external boundary and 
the plate. The closer the external boundary layer is to the plate 
the greater the influence of the intersection of reflected 
diffuse lines with the wave on its intensity at the boundary (the 
intersection causes bending of the wave and thus reinforces it). 

Since ^ calculation of intensity of the incident wave for P e /P Q 
is anyway iterative, that phenomenon is not important during 
the initial stage of evaluating shock intensity and Its In- 
fluence is anyway expressed by the solution of the field flow 
during continuation of the iterative computation. 

Just the same, it is not desirable that most of the re- 
flected diffuse lines reach the external boundary upstream of 
the incident wave, since the influence of the wave’s bending in 
the flow field will then not be properly accounted for; on the 



(C-7) 
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other hand, it is not desirable to have the external boundary 
too far from the plate, so as not to increase the number of 
gridpoints unduly. Experience has shown that as long as no 
more than 25 $ of the reflected diffuse flow lines intersect 
the external boundary before they meet up with the wave it- 
self, the solution is not much affected. Based on this 
criterion, the width of the field, ly, is determined in 
appendix D * . 

In addition to the factors mentioned, the influence of Ax 
on the accuracy of the reflected diffuse flow lines must be 
noted. Ax must be at least one order of magnitude smaller 
than the length of the reflected diffuse lines zone so that 
their influence be a close approximation of the real situation 
(see appendix D'). 

During p ractic al computation an initial evaluation is made 
of = /p /p Q , which is corrected in accordance with the 

results obtained until the desired value for p /p is reached. 

o 
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Appendix ,D' - Optimum Evaluation of Cbmputation Grid Dimensions 
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The computation grid must include the region in which the 
majority of physical changes in the flow field occur. The grid 
density must be sufficiently high to permit observation of the 
main phenomena in the flow field, but it is also desirable to 
keep the number of grid points as small as possible so as not to 
place too heavy a load on the numerical computation. Also de- 
sirable is that the ratio of of the dimensions Ax/Ay not be too 
large for reasons of computation stability. 

These considerations will be brought to bear in the evalua- 
tion of Ax, Ay, lx, ly. 

a. Evaluation of lx 


Based on experimental results and on theoretical evaluations 
(see also paragraph 2.1.2), it was found that the region in which 
the principal changes in the flow field occur is about 1506 o 
long and is divided into two more or less equal parts on both 
sides of the point where the shock wave enters the boundary layer. 

Therefore lx = 1506o (D-l) 


During earlier computations (at the start of this study) it 
was found that, based on this evaluation, the location of the 
boundary downstream is far enough from the zone of interaction 
to fully justify the definition of a boundary condition of zero 
gradients in the x-direction. Also investigated and verified 
was the existence of the adjustment condition 
(in addition to at that location). 


Op 

^ r “> 



82 



From the point of view of the region where changes occur in 
the flow field, ly must be or the order of magnitude of 2 - 3<$o 
(as obtained from previous theoretical and experimental results). 

We do have to consider an additional factor in this study, 
though, and that is the proximity of the incident wave to the 
diffuse reflected lines, which stems from the closeness of the 
external boundary to the plate; for that reason we will employ 
the sketch. The shock wave enters the upper boundary at point /88 

x-y at the angle £ and its continuation reaches the plate at 
point x (with approximate assumption that the wave will not bend 
and will arrive at the plate). 


The separated region stretches between points x sl and x^ 

and its length is Ax . From the end of the boundary layer above 
this section emerge s the diffuse reflected flow lines (compression 
lines at the ends and expansion lines in the middle). 


Based on previous theoretical and experimental results, it is 

known that the largest part of the separated region is located 

before point X g (because of the influence to the rear, exerted 

by the pressure rise through the boundary layer); for the purpose 

of an approximate calculation we will then assume that 2/3 of 

the Ax section is located before point x . 
s s 

Since the density of the diffuse flow lines is much lower 
at the periphery that in the center of the reflection zone, we 
will only concern ourselves with the central 2/3 of the re- 
flection zone in this evaluation. 


Let us also assume that the angle of reflection of the 
diffuse lines is equal to 0 (actually the value is very close 
to e, because of the small local flow angle which is on the 
order of just a few degrees). 

Based on the foregoing data and evaluations we will now 
calculate the distance of points x. and x-.. from the front end 

1 KX 

of the field. It must be ascertained that A is positive or, 
at the least, of very small negative value so that all, or the 
majority, of the reflected flow lines will go to the external 
boundary downstream of the shock wave. 


Based on the foregoing evaluations 


O.S i 


(D-2) 
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- X 5 ,r AXs 


we then define 


_ 7 , 


■ £ “ Xft » ” 


(D-3) 


and get 


X = £o. £ i-* “ o, 6^} ( o, i> V^Ks)J + 
“ 0. S'iy - o.^T AX S 4 


l 

V 


1 was already determined, previously as 1 = 150. 6 o 

X X 

we designate xc/jv x^/r^xt. 

AXi/<fo = X h/s>* ^ 


(D-4) 


and ultimately obtain = 5- 


x A - ^ 


(0.-5) 


Prom computational results for 2< Mo< 5 and 1 < P e /P Q <3 /89 

(which agree with 10° < I <40°) the following orders of magnitude 
are known for the separated zone: 

Ax =(20 4 50) 6o (D-6) 

s 

Now we will investigate the influence of ly on A in the 
region (1 - 7)<5o : for this purpose x 1 and are calculated 

for 16 different combinations of 6 and 


Ax for which 
s 

Ax s = 206o; 306o; 406o; 506o 
0 = 10°; 20°; 30°; 40°; 

and the results are tested for the following values of 1 /<So: 
1, 3) 5, 7- The results are sketched on the next page. y 
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X % of cases that 

'fulfils the indicated 
region of 


£ >o 


0 < ~ < -o.l 


0 - < -o.X 

£* A£ 


O < — <- U. 1 
^* 4 . 


, , r>o<snT , ts from 64 combinations 

o^l S for which the value A was calculated 

according' to formulas D-5 


Prom the computational results it appears that in the 
region 1.5<5o <1 < 3<$o the value f increased greatly for a 

y 

A that was either positive or slightly negative, while in the 
range 3<5o <1 < 7<$o this increase occurs much slower and more 

y 

gradually until it reaches an asymptotic value. The results can 
be summed up in the following table: 



r'N 



1 

3 

5 

7 

A > 0 

10% 

50% 

68% 

81% 


30% 

89% 

95% 

100% 


It appears then that l y = 35 o supplies effectively the re- 
quirement that no more than a quarter of the reflected diffuse 
lines reach the external boundary before the shock wave. The 
improvement available from 1 = 560 or 1 = 7 < 5 o carries little 

y y 

weight in consideration of the needed increase in gridpoints 
( 1.7 times and 2.4 times, respectively). The optimum choice 
will thus be - 

1 = 3$o (D-7) 

<J 

c. Evaluation of Ay 

To be able to detect changes along the width of the boundary 
layer Ay must be an order of magnitude smaller than 60 ; so as 
not to exceed this principle the maximum value of Ay will be 

Ay = I/860 (D-8) 

Note: Because of the increase of the boundary layer by Ay < O.I60 

in the interaction zone proper. 

d. Evaluation of x (For Uniform Grid) 


There are two requirements that limit Ax and they are: 

1 ) . For the sake of stability it is desirable to maintain 
a dimensional ratio of Ax/ Ay < 20 and because it has been 
determined that Ay = 1 / 86 o it is required that Ax < 2 . 5 <$o 

2 ) . Ax must be an order of magnitude smaller than the 
length of the reflected flow lines zone so that the lines may be 
detected. That length is, according to its evaluation in para- 
graph b., about 2/3 x (20 - 50) 60 = (13 - 33) 60 

and therefore 1 < Ax/60 < 3 must be valid. 

From those two requirements it was determined that 

Ax = 2 . 6 o (D- 9 ) 
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e. 


Evaluation' of x (For a Non-Uniform Grid) 




For a non-uniform grid it is required that in its center 
Ax = 0(Ay) be valid. In this study the value chosen is Ax 
min = 1.6Ay = 0.26o (D-10) 

According to the ratio of 1.1 for growth of grid meshes 
towards the periphery (see considerations in appendix G') we 
find that Ax < 26o remains valid in the zone of reflected 
flow lines. 
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With the coefficients a 3 a , b, . b 0 , b the choice of subjects 
increases x y x y 



The following relations exist (from now on all length values 
are dimensional): 

According to appendix A’ 



(E-4) 
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also valid is 


X ! (E-5) 


£-X 3 ** So 


(E- 6 ) 


Those are the three equations with the unknowns 1 , x and 
60 and the procedure for their solution is as follows : x u 

jt« y, I J J ' 

After extraction of <5o and x we get - K, V - — 

o ft./. * ■ 2 - 


for 1 


We designate and get a quadratic equation 


•sXVKix - -2 f< s 


(E-7) 


From the two solutions we choose the practical physical 
value and get 

^ a »K - { ^ K(j J+i6xs/k - -<) (e- 8 ) 


while and x Q are determined through (E-5) and (E- 6 ). Through 
(E-l), (E-2), E-3) and <$ 0 we determine 1 and also Ax and Ay • 

tj 

(after substituting the coefficients in those formulas). 

The number of gridpoints in the y-direction in both grids 
is (E - 9) 

The number of gridpoints in the x-direction of the uniform 
grid is . f . (E-10 ) 

The number of gridpoints in the x-direction, in a non- 
uniform grid, is based on the mesh length in the center of the 


field being Ax 


The mesh length increases from the center 


of the field towards its periphery according to the formula for 
a geometrical series with a growth ratio g. 


So we get 

when n is determined by 


1*- - t 


M, = ( E-ll ) 

(E-12 ) 


9 _ CX jl 


Jc&ii cj - I ^ 

and its solution is 


A jCuiv 


Tt- 


l 


(E-13) 
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p 


According to the previously determined considerations the 
following coefficient values were chosen: 

1; (e-14) 

<V J J 


and then we get the number of gridpoints ' 

ri s -r[y 

The number of points was limited to this value (76 x 25 = 
1900 ) by considerations for an optimum between high accuracy and 
reasonable computer time (one computer run without calculations 
of coefficient stability for 1900 gridpoints takes 1.0 - 1.5 
minutes and at least several hundred runs (iterations) are re- 
quired until convergence is reached). 
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Appendix F' - Computation of Flow Characteristics in the Shock 
Wave Direction at the External Boundary * 
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Given a shock wave with intensity Sh = P p /P E 

Shock wave 
External 

It is required to calculate the boundary - 

flow at F as function of the 
characteristics at E. 

From the Rankine-Hugoniot relations we get the angle at 
which the wave enters 



’ « i/:-} 1 yf-8 




— «.? r 7 / < \ • ' V *' !>‘rl 1 

i -~t* Jl 

Velocity components at E: M?, s U x zs 

across the shock wave .*/• ,, 

V- £i .•? bib- &•*</* -»4' sui 

parallel to the shock wave 

, S . 7i 1 

(; r-;y.v.g 1 


(F-l) 

(P-2) 


(F-3) 


across the shock wave 
parallel to the shock wave 


V t -- « V- 


(F-4 ) 


Velocity components at F according to the original directions 
(of the first axes) 


in the x-direction tip s 2,w-</r + 

in the y-direction Vf ■VJ.caS; - v$, s SJU3J 

We also get 


r , t. f l£ 


a( ar-y 

and from the state equation we get 




J** £: 


r fv 


-S TV 


(F-5) 


(F-6) 

(F-7) 
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Appendix G' - Formulas for Finite Difference Computation 
a. Uniform grid, second order accuracy in the x and y directions 

According to Taylor’s development around m_ 

~ /v- Jy f | 

m- 1 i* - ! m+i ! 








L 






m-t 


5 . 

t 

1! 

L - k t ■■ 

oi'* 

& 

h l 
a i 

r: 

a*} I ^ 

f 

.C 4-K ^ 4 

j Wl j <<t T 

M- 

,3l 

r; 

* i> 

£f 

f« 

+ AV 

o<:V T>_ 

subtraction 

we get jv v 

* 


h " f 
«ilv 

1 

■ 4- o ( k * ) 


by addition we get 


(G-2) 


U, Y auUJ.bJ.Uil , (■ T f It 

^0(^1 

For a grid with uniform spacing in the x and y directions 
we get i , 




^ 

*7 




/' 

" <1 f- { > ” f~ W*-i , v». / 1 v •. 


'EX j W\,„ 

°4- 

'j 


A l •.* 

W I'.-W - f M , -> 


*( C; 

* ( ^/ v j 


1 ( G— 3 ) 


n> K ~b _ “«2 / 1 t\ 

— , - — — — 1 « — tOI W* j 

' ok 4 Ly K ' 


m * 

<&_V 

■***' 


-P *n , r. ^ rf u\, * - 1 


(j~*V4l t l».4l 2 ~f fr-.4l.rt-l) “ (fm.w-U " 'f 

'VAx k.y 


+c (^) 


b. Uniform grid, fourth order accuracy in the x-direction (and 
second order accuracy in the y-dlrection) ~ 
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According to Taylor's development around m 

lx U. U. 4 tv . / 

M-V V-\ Vw WAV UA4U /(G-4a) 

-M v -9 s 

k|-4 £|V 4 + o(Vj <a-«b) 

tKM . f.* W« 4 ♦ £/!. - ^ *« + ; <0-40 

|(a-»d) 

41 . . -U - 2^zli -■-> “* - 4-u ”-* -1^ 4-0 CM Y 

^ 4A - W V (G _ 5 ) 

f _ - fm-2- -v^b-f^-v -so-fbA *-H»4w +i -f*A_^ 4- O (. k q ) 

T <fa_V* v 

For a grid with uniform spacing in x and y directions, with 
fourth order differences in the x-direction and second orde 
differences in the y-direction, we get 


(G-5) 
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(G-6) 



= 5 V f* -*■»*+! ~ * ? j"K*AI,A4» ~ f ^ + i.,M4l) + 

- (fu AM - 2 j-M -»,«-» + J flwv 41, *H “fv*4 j 4 0 »Sy 


n- 


Note : when h is one order of magnitude smaller than h (for the 

y h ^ 

uniform grid of the computation' ^ = 16), a higher accuracy is 

h 

obtained in (G-3) for de.rivatives y of y, which is not expressed 
in the finite difference equations because of the lower accuracy 
for x-derivatives . 

In contrast to it, the accuracy of x-direction derivatives 
improves in (G-6) and reaches the same order as the derivatives 
in the y-direction, so that in general the accuracy of the 
difference equations improved. This improvement in accuracy 
is particularly convenient near where the shock wave impinges on 
the boundary layer, where relatively strong gradients of the 
variables exist. 

c. Non-Uniform grid in the x-direction, second order accuracy in 
the x and y directions 

According to Taylor's development around m 


Ivv* 






■^4" 4 

b 1 ** 

\j * ^ 


" 



-Li’* k ** 

* " f K, + — 

1 oi 

4* ■ 
T M 

4 * ( 1* s ) 

(G- 

-7a) 


1 — rpk* 


ct * t\ * , 

1 y & •*» 

* -f m f b ( u 9 

I (G_ 

-7b) 

Through 

1 * (bj we 

get 







f \ s 

4-U -.1 

*" 4 •*•-') * °i . _ ^ 

J(G- 

-8) 

Through r . 

t ** *7 

■lx fa] we 

get 


rrr 

li 




f'* r . 

*4 4 

1 ~ «~P w>) ~ ( ~ i£ U* - 1 \ 

I 


when 

the errors 

T 

are ■' 





J 


fM e 

<jUU*<\) 1 

f- T 

W 


sp>o(h i ~ c {) |; 

(G- 

-9) 

K 


[■PiliVA 

■ 






p 

For q = 1 the errors in the two derivatives are 0(h );^to 
preserve the order of magnitude it is .necessary that 1 - q = h. 
In practice, the accuracy of the derivatives becomes impaired 
when q <1.1 and so we will choose q = 1.1 for the difference 
computations of a non-uniform grid. 
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For a grid with changing spacings in the x-directions 
unifo™ onef in the y-direction, when the differences are 
computed to a second order of accuracy in both direction , 

get 


and 

we 
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Appendix H* - Detailed Finite Difference Equations 
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rs 


a. The Brailovskaya scheme (for uniform and non-uniform grids) 

The complete difference equations will be formulated 
further on, in their general form and to a second order accuracy 
(without artificial viscosity terms). Let us define the symbols 
as follows 


1 

J 

Kon^ni^orm grid 

Uniform. Grid 

1 

J 

A X Wi /A X W»-8 

4.0 

II 

£ ( AX tAX^.i] 

AX 


A< 3 



The general form of the equations 

D W 

First integration stage 




= S 


(H-l) 


vv *«* = * (fJL* - 

^ + At- S w, >v . 


* 2 . 


Second integration stage 


(H-2 ) 


- aJi 1 

f/p 4 - - 



if h *' - f V Hi 

c2 

[\ *44, H 

W»,vt J 


\ M.V4 K\*l |« / « J 

- At 

r & lT ‘ - 

r tT ' 

j- 

At ■ 

Vn*W 


The equations for fourth order accuracy are written in 
similar form according to appendix G ' paragraph b. 
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The expressions for the equations of a uniform grid, to 
second order accuracy, will be detailed further on in the form 
in which they appear in the original Brailovskaya [8] scheme 
and also in the computer program. (parallel terms for a non- 
uniform grid and for fourth order differences, were developed 
for the computer program in similar fashion, by means of the 
definitions given in appendix G’). 


W e 

«V* 


'/ ' 



fix 


Wz 


j 

v/± 

E 


Wu 


(H-3) 
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rr 


Vv'l 

• . 

p*j3 W v 


|° 4- V/ 2. • U 


s w/ 


WiV 


u(6^p) 

_ 


[/.{v/il -4- p) 

Im.vi 

’ ' 


w'i 


j>w 


Wj>; W 

* 

pj- 




Pj 





(H-U) 


(H-5) 


c C J_ J. 

Kv'b X 


TVy 




v'Tv 


f*r d>X 


7 






Zxy 


(H-6) 


I? * UT «1 ^ 
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r \ 


W s + *; t s + 7 r s t 
hrst-f r s + rirs f r , s 
o 


-Z = (H-6) 


When 


(H-7 ) 


/100 


S * xs 

^• 3l ' *’ *}|;l^)‘ ’diJ,[ - >**»*'.- (v/-h, w *, - vc., lW .,)J 

• * S3 ^' • 3 ' t ^^ A ^)" Jiv »{( h .--• Xi / *-.- i - c .«- 5 j 

5 i3 u ^.)] 

5 3^5 r - - Ll k ^ A - _ J * f ( / 

J - UjV ' ow ;- JJ ^/ A - srt *.( U WS 4 , | 1 ) , 4 i - W ».. i,«y : ^ wl | ». i ( u Wk *, lW , l - U.^^.JJ 

^ ^ •**) 'aJvT 7 [(^. *j A **- v ')( Twa ' «** -,.1, .) [x, -T„ . , ^ jj 

5 I,-;* z fti./u u W . 3 TV .. 

i^lT '««; ~ jt ^ c ( Ip .-..,* .,,„) 
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( H— 7 ) 


^ ' 3 - ^.-^)] ' 

= jiW^ [(^‘-. , »*» w i-' .v*A* U.J 

S i 

^ ^ r ' /TTvTtT^ [^’*' - M ,% * ' / ^» ». *• (^<«»> ( WM- l/u -J a **.,,* \A*-(,» (u ^»,v**C-Uh*.|,M.i)"| j 

^ A'iy (^ AVA ‘a , j)' sfXlnWy ^ < M w '<l,n'^i«*t lk , ( V'U4 J| M«I - l/i*»4|,M -|) **^A twJ'pU M.l,^(l/| lt | i i,(| - K»-t,n-*Jj 

J i 

~'5^(| A ‘ U 'i3c} “ [^ •».*»*• Ul ". IMI ( VC«|,W-U -^>M,lW U Km-« (\/v«4l,v*-|- 


b. MacCormack scheme (onl y fo r a uniform grid) 
The general form of the equations 


iw + a_£ ^ "yr _ s _ \£j + \s± 

1-i T>X TiM " ■ -OX. T>^- 

and in another form - " — 

! Ti£ 7 ^ q . . 

j ”&i T>X T>M ~ 0 

! I 

6 = &-il 

Here follows the general description of the scheme: 
The first integration stage 


(H-8) 


(H-9) 


i +■ 

V/ = vVi 
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when 


hAi =■ Va 4- vw A 
- vs * v\ A 
wvA - 

V\ A r Vwo J ( ^ - Vw<J y ^_ J 


(H-ll ) 


The second integration stage 

t-M t“» 


FT. 


, T + U, .4 , fc-M X 4 , t-*| \ 

w v* t * r a l w ^ + w '~~ - 

| Wl - Cw + vw A 


(H-12 ) 


when 


(H-13) 


j I'j T V\-* v\A 

V^ol - VAO A ^ Wy A 4- | j 

v\A =■ w . 6il ( vn A A *, 


Note : the function mod is defined by 

Wod^,^ x - tA/T ( 2L) , ^ (h-14) 


W »f 


(H-15) 


Ul.K' 


' f ' 


WA 

J=M 


Wx 

f V 



E 

k.Vt 

wU 


f) 




■ r - 



o 

P+f w v 


4 

c T* 

J.UV/ 


R« 

Tyv| 

u(E+f) 



NK 'b T 

. fN ^ + «V ^ 


(H-16) 


(VIlA 
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(H-16 



— 

w ^ 


I 5 4-1 V^U. 

JA *v» ,** 

Wav' 


ulwH^) 





o 

A [ax ( * U '*-•!*) * 

jjnu^»,»4 -T w ., ( J-t- — u {y**,*" 

* ^■'•''j^iy ~ ( ^*>* 1 x'<] +j £ ~jl ( W*m,«. ~ l^*i i») J 


'r 


0 

j> u */ 


Xoty 


f?4t 


v{e*t) 

kvi.n 

. 5 pT^' * ur ^ y 


(H-17) 


vVi 
Wi*H- 
PMViV 
v(vjh4 »*) 


— 


_£■ [^(v^,^ 1 ~ -JT* W***-*i|i*'’^‘*'' , t l *)*j 

p, ^ [^I- ^ ♦!“ ^.n^J -4 ^ in . V» I jXvj C 1 ^ 1-^ , ^* *• “ Uwi.«-«) *^2x1$ 


During the actual integration the computation is split 
up into x and y directions. For all KI runs in the y-direction 
a computation for the x-direction is made where KI = INT 
(DTX/DTY ) 


DT/-o.r Ax/jr* j>Tr= o.r^y/ [cxo(uo.irH,)] 


(H-18) 
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A set of computer runs for the y-direction is I while for 
the x-direction it is II when 

J4 - jj/T ( I / Ki) (h-19) 

Details of the computation procedure 

1). Computation in the x-direction - 

First the values A, B are calculated according to the defini- 
tions and then integration of the difference scheme is carried 
out per time period ATx . 




(H-20) 


when. _ 

(vU*~ V*-* I 4 - Xa/T 


"J 


("A) ** 

* r \ lA ■+ m - ivu -Xa/tC 27 

I T fl + St/T] - — ; 


f/f- tVI .Vtvl 




-} 


V (H-21) 


This arrangement generates alternately forward and rear 
differences in the x and y directions. 


From the results of the W values the values p, y , T, p, V, U, 
are calculated and from them the values of A. 


Finally, the average derivative in the x-direction is 
calculated, for use in the KI runs to come in which only inte- 
gration in the y-direction is carried out. 


when 



(H-22) 


(H-23) 


Note: Experience has shown the desirability of the definition 

[F ] = 0 for the first run (because the calculation of the 
m,n 

derivative is not accurate) and afterwards the computation stage 
for the y-direction is performed. 


102 



2). Computation in the y-direction - 


/10H 




First the values for B are calculated, according to the 
definitions and then the stages of integration are performed, 
according to the time period A TY. 


First integration 


when 


hj - <*4 1 -Xa/t 


(H-24) 

(H-25) 


This arrangement generates alternately forward and rear 
differences. 


From the results of V.’ values the values p, p, T, P, V, U 
are calculated alnd from them those for B. 


Second integration 

j'j * -+/ij 4 -X WT ( *2- 

= h+X - £a/t( l/2.)*x 


(H-26) 


when 


} 


(H-27) 


This arrangement generates alternately forward and rear 
differences, which are opposite in direction to those of the 
first integration. 

3). Final calculation of the variables U, V, p, T, p, p 
is performed from the final values obtained for W. 
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Appendix I* - Details of boundary conditions and other computa- 
tions in finite differences 


a. Boundary Conditions at the Entry Section 

According to formula (3-19) the relation f ^ (1-1) 

is valid and in that way the suitable n is obtained for 
all gxidpoints in the y-direction. Prom equations (3-16) to 
(3-18) we can obtain a , -p » ^ - as function of 






+ *l» (i_2) 


T "'" ' & tr »•{*- O] 

J> 4 tj _ ([ 

4-h <T r ~ Ltt,n) 

VVjm - 0 


(1-3) 

(I-M 

(1-5) 


and from the equations of state and of viscosity we also get 


b. Boundary Conditions Along the Wall 



I 


l\v> 


Velocity components 


from the non-slip condition 

u*,. 


( 1 - 6 ) 


when for suction 

and for injection Vv* & 

Temperature must fulfill the adiabatic condition 




and we will assume the following relation (which fulfills that 
condition as : /Ti*,*, * Cv+^vit. (1-7) 

The coefficients a and bvwe will find by substituting the 
values//-^ ^ 4 v \ into the equation 


« 0*t\) 

• (Sv) 




-ft 






T WN, V + *■ ’ 
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From the solution of the two equations we get 


but 'M|» 


— CX. 

V s - 




wi,l * Tw>,} 


J> A 


\ A r 


- ti-r --It- ' 

' 3 Tw ii- j> r ‘"'ij 


and therefore 


«T _ ^ T _ X T 

! *W\,| ~ 3 *Vn iV “ 3 (1-8) 
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According to the considerations in paragraph 3. 3 .1.2, 
pressure will be computed for the condition » | 

<&Jl1 ~o 

Similar to the temperature computation we then get 

p> _ fi p _ JL p 

r *\,i " 3 rw i3 (1-9) 

and from the equations of state and viscosity we also get 

Here we will bring up a few more methods for computing 
which are basically more accurate (without the approximate 

assumption of _5£ _ q) resu ;Lt an unsteady solution for 
numerical computation: 


1). Linear extrapolation - 
Assuming the relation 






we get "* 1 

2) . Square extrapolation - 

Assuming the relation j f'^x* + + C. 

we get 

3) . Development of continuity equation 


(1-10) 


(1-11) 


H vuJ-A +- I'H » o (I _ 12) 

xb ^>X •'5^ \ S K V n ^ 

* u i ^ 

T.T -1 1 * 3 . _ • « 


We develop the expression around point; 
while neglecting the second order 



[ 

><*•) 




L_ 

L 


/~N 


and finally we get 

(t) (i-i) 

f ^i> ~ pH, I 






A.t 


+Vi 2. ~ ^ 

J«W L - 


of 


(t) 


(fc) 




(1-13) 
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(t-o 


p w = 

J **• 


it 


4 + A t 


L o 


vn 

Hn - 1 1 1* 


v*' 1 
+ ^|V 


^4 ) . Development of the approximate equation of continuity - 

This is the approximate formulation of the continuity 
equation _ • 

U i= - 

-fc-t I (1-14) 

and after development of the expression around point { \ 
we get V 1 ' 

(i.) n .\ tk-0 

>. _ p*,x • W (1-15) 

- f- zq; 

c. Boundary Conditions at the Exit Section 

Based on the gradients of the variables in the x-direction 
being zero, we get - - ■ 


VV,,k * " 

A - l*M-1 . * 

Tm.w -Tm- 1 ." 


( X— 16 ) 


d. Boundary Conditions for the External Flow 

Computation of the characteristics for the external 
boundary: * 


// r— 

7 

/ 

/L 

/ 

/ 

/ M 

\tl 

V 

/ 


’/v/-1 

The characterist: 

if i 

H * 

Lcs are c 

r — i 

a ^ 

omputed 

i i 
l4« 

from 

the 

next to the 


last line vs. the last line and their direction 3m is defined 
by t 

*9-,,., ) * Yfeir) < M » 


The value of a random variable F „ will thus be determined 

m,N 

as a function of its values from line N-l 


^jKI 


( A * u ) * ~ l) 

C^i.X-AXl «*• <iXu) 


(1-18) 


/1Q7 
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when the lines of the characteristics bisect the vertical grid 
lines, the definition of F M will change accordingly. 

The values at the external boundary, at the first point 
downstream of the shock wave, are defined according to equa- 
tions (3-26) to (3-29) in paragraph 3. 3. 1.4, so that the values 
of point (m ,N) ahead of the shock wave define the values of 
point (m + f , N ) behind it. 

S 

For reasons of computation stability we compared the values /108 

of point (m +2, N) with those of (M + 1,N): there appears to be 
s s 

almost no effect on the results but it does tend to prevent 
fluctuations in the solution. 


The value of Vm +1,N is not computed from the transition 

equation of the shock wave like the rest of the variables, be- 
cause the changes in this variable are of the same amount as 
its actual value and' such a computation would cause severe 
oscillations in the solution. This value is therefore computed 
through extrapolation of V from the field by means of the charac- 
teristics method, similar to calculation of the rest of the points 
at the external boundary. 


e. Boundary Conditions on Both Sides of the Shock Wave in the 
Field, for "Non-Continuous" Computation. 

This boundary condition is computed according to the 
following stages: 

1). The first iteration 

gives the point of 

entry into the field 

(M .N ) and its 
o o 

intensity Sho from 
which the entrance 
angle is calculated: 



tv 

4 






i(ik 


. 2T-l\ 4*1 *f 

Mu 



(1-19) 


2). Computation of the initial location of the shock wave 
in the field (assuming a straight shock wave) 


h) - Va <. ( h. - h-j y C'J.) ( I_2 0 ) 




3). Setting of the upper boundary of the shock wave Mu(N) 
and of its downstream boundary I , , . 

I h) " 
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4). Location of the lower boundary at n so that it is 
close to the sonic line where M = 1 (through s computat ional 
test ) . 


5). Computation of values upstream of the shock wave (for 
points M (n)) through extrapolation of values from the upstream 
flow. u 


6). Computation of the values downstream of the shock 
wave (along points from values calculated for M u (n) 

from formulas for the transition of a slanted shock wave, which 
are given in paragraph 3. 3.1.4). 


The value for the velocity component V downstream of the 
shock wave is not computed that way (for the same reasons as 
discussed in previously in section d.) but through extrapola- 
tion from nearby points at the downstream flow. 


7). After correcting the computation for the entire flow 
field by means of the difference scheme, we find the new local 
gradient of the shock wave from the local pressure ratio. 

on both sides of the shock wave boundaries from the 


formula 




= SK J [£1 


t?-\ \ >-»r 

/ a .y 


M 


} 


( 1 - 21 ) 


8). From the local gradient the location of the shock 
wave is brought up to date 


J ^ (1-22) 

. 


afterwards one continues according to paragraph (3) to (8) and 
so on repeatedly. 


f. Additional Computations (Based on the Definitions in 
Appendix A* ) " ’ 


1). Mach number 



V >-l) 


I 


(1-23) 


2). Flow function 


(1-24) 



‘/'ki-i.h-i 


-l- 





1 & 

jjjT C < 11 "* 

f ^ VV\4, Hm) 

AX 

^ *2. ) 
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3). Friction Coefficient along the plate - 
Per definition 


(1-25) 




We will calculate v. , in the following manner: 
we assume that the following relation exists next to the plate 
by substituting the points - |va*\ we get 

V-K.= «•» , (*,v) , l*'') 

u* 4 a- 0 '-o' +0 ^'^ * 

from the solution of the three equations 

•-J Ay ' 

and finally we get , , 

C ^»n r 'RcTT'Sm / Vw *‘ + Mm 'V (i_26) 


4). Heat transfer coefficient along the plate - 

. . <w v . » . tA ~r\ 

Per definition 


we calculate 




c ^wT 

»» » » ’ 
P, 

like 


. ■ / 

+X», v) 





and get 


Finally, we get 


(1-27) 


( 1 - 28 ) 
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Appendix J* - Comparison of Stability Criteria for Computation 
of Viscous Plow According to Various Sources 


Based on various approximations in the development of 
linear equations for compressive and viscuous flow, different 
At criteria were obtained for assurance of stability in the 
computation. Some of the most important results (in dimension- 
less form) are shown below: 


a. According to Brailovskaya[8] ^ ^ 

without the energy equation v *£ 

b. According to Brailovskaya[8] with the energy equation 


Pr £<, Kj'Cy-O _ 




A-t ^ g 

c. According to Skoglund, Gay[9^] 

At k — ^ e _ 

d. According to Carter [13] 

A t £ - fr 


o.s f, g P 




(J-l) 

l 

(J-2) 

(J-3) 

(J-4) 


If we set the non-dimensional values (the approximate 
values) | Pr ~ 0 . ? rw<ll 


we get 


o.ii-S 

0.03,5 M„^ ( 
- <3. 1 % Ke. 

f= c - O t XS e.e. 


(J-5) 


For 2< Mo< U all the expressions offer approximate values 
in the range F/Re = 0.1 - 0.2. 

In this research only the criterion of F q was taken into 

account because it is based on a minimum of approximations of 
the basic equations. 



Appendix K* - Development of Coefficients of the Differential 
Equation Remainders 

Arrangement of the general computation with difference 
schemes for the t-th iteration 


For the Brailovskaya scheme 



derivative 

vA • w 1- ' - At 

nnrxa J 

i i . 

P derivative 

w = W* -At 

I imuj 

" J 

: scheme 



/ - . , r derivative -t-* 

iv/* = W 4 ' - Ai ... iUj 


l 


derivative - 


w 


“ n 11 Til J 


(K-l) 


( K-2 ) 


when all derivatives in each scheme are calculated as described 
in appendix H*. The remainder of the differential equation is 
defined by derivative 

lt-« 




(K-3) 


According to the rear differentials we get 

Vt-» 


(-busy-' _ y/t-'L 

J At 


From the previous definitions 
(for both schemes) 


(K-l!) 

(K-5) 


derivative 


] 




At 


In the same way of making use of the approximation of the rear 
derivative of 6 w r- can be exchanged for TT t' 

fit W t , W * 


Then we get 


£ 


*-• 


At 


(K-6) 


This expression provides, of course, only an order of 
magnitude estimate for the remainder and not its exact. value, 
because of the approximations in the computation procedure. 
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Appendix L* - General Evaluation of the Computation Error 

We will perform an evaluation for the following reference 

conditions: p /p = 1.4; Re = 3 x 10 5 ; Mo = 2; 

*e o * x 

s -2 

Reference dimensions are: x „=x =5x10 

rei s 

<5 = 4 x 10 
o 
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Prom this we get: 


4- - 5x1 o' 

L'* .i "" i > * n 


r 2 


^ -- <fv/ x 


7 00 

-4 


=7x10 


-5 


Ai'o 

« - bx . .-2 


8xl °-o %6 cl O'* 2 7 




M 




cT > /x j- 

5x1 o“ 5 

' -p 

5x10 * 

0, £ A*t CFL 

5xl0~ 8 


’ 7xl0" 5 

^51a/ 

t v ^ 

- 1 



= 1x10** 3 = 


7x10 


-4 


e( io” 2 ) 

< 10 " 3 ) 

o ( lO** 3 ) 4-o(lo’ 4 ) 


The truncation error for the Brailovskaya scheme is 
determined through E T = o(At, Ax 2 , Ay 2 ) 

If we take into consideration the influence of the 
oscillatory solution fcomponent as well (because there is no 
significance to the dependence on At in the steady solution), 
we shall see that the truncation error (which is of the same 
order of magnitude as the reminders of the differential equations) 
is E t = (10- 3 ) - (10- 4 ). The order of magnitude of the variables 
themselves, however, is At times larger, in other words 
( 10 - G ) - ( 10 - 7 ). 

This analysis does not take into account the value of the 
double derivatives At, Ax, Ay in the complete formula for trunca- 
tion error which, in certain regions, may increase to several 
times that amount. 
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It is also known that near the boundaries, and particularly 
in the region where the shock wave enters, the truncation error 
is greater since the local deviation between the accurate flow 
equations and the difference equations is larger. 

As we will see from the results (see Pigs. 10 -Ik and also 
appendix M* ) an average order of magnitude of 0(10- 3 ) - 0(10- 4 ) 
is obtained for the errors of variables which indicates that the 
lastmentioned factors exert significant influence. 



Appendix M* - Evaluation of the Computational Error (Through 
"Comparison of Kesuits from Grids of Different Density ) 

A comparison of results was made between those from the 
regular grid (run 001) and those from the fine grid (run 041) 
for the following parameters: 

x s = 0.05m; P e /P Q = 1.5; Re x = 3xl0 5 ; M0 = 2; 

s ^ 

Regular grid parameters - Ax = 7.97 x 10 ; Ay = ij.98 x 10“ 5 

Regular grid parameters - Ax = 4.10xl0 _i| ; Ay = 2.51 x 10 -5 

(The exact ratio for the mesh of the grids is 1:1.944 in the 
x-direction and 1:1.984 in the y-direction) . 

The results obtained for two characteristic points (for 
which the location comes closest to being alike for both grids) 



Type of grid 

Non-dimensional values 

Location in 
the field 

U 

f> 

T 

M 

A/ 

Before 

the 

shock wave 

Regular grid 

0.1377 

0.1730 

1 .0250 

19 

2 

Fine grid 

0.1493 
0.1481 * 

0.1751 
0.1751 * 

1 .0293 
1 .0294* 

36 

3 

After 

the 

shock wave 

Regular grid 

’ 

0.2004 

0.2585 

1.1171 

55 

2 

Fine grid 

°.1973 

0.1957* 

0.2522 
O.2522 ■* 

1 . 1 209 
1.1211* 

106 

3 


*Since the points from the regular grid and from the fine 
grid shown here do not match up completely, an interpolation was 
carried out from the fine grid results to the appropriate points 
in the regular grid. 

Let us assume that the computational grid errors relate to 
each other as the square of the grid spacings (and that is 
based on a truncation error of the order 0(h 2 ) of the difference 
scheme ) 


n 
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(M-l ) 


r-\ 


f c - for the regular grid $ c - f 


fp - for the fine grid 


- i ( h A) 

Prom here we get the exact value of the solution 

! ^ _ 6 ~f /■= ~ -p*. 


= 4 


/Ilk 


(M-2) 


That way we will find the exact values for the points 
examined and any divergence from them in the two grids. 


1 

I 

1 

1 

Type of grid 

Non-dimensional values 

Location in 
tha field: 

U 

' 

P 

T 

H hi 

Before 

| 

the 

shock wave 

Regular grid 
Fine grid 

0.1510 
= .0139 
— • 0035 

0.1758 
=.0028 
-afe =.0007 

1.0310 

=.0064 

^-r P =.ooi6 

19 2 

35 3 

After 
; ^"the 
jshock wave 

Regular grid 
Fine grid 

— . 

0.1941 

=.0063 

<=AJ_ = .0016 
* 

0.2501 
^ = » 0084 

= o 0021 

1.1220 
^C & = . 0049 
tX? =.0012 

55 :2 

106 3 


It turns out that the errors in the regular grid are on 
the order of 0(10- 3 ) approximately (in the fine grid they are 
smaller by a factor of 4) . 


n 
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Appen dix N 1 - Comparison of Errors and Computation Times in 
Various Studies ~ — 


A table of comparison was prepared, based on various 
studies made previously: 
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Paramete 

Compared 

Lax-WendrofT M..ebhod 

' - J 

Brailovskaya Method 


.MdcCormack 

[67] 

MC 

Skoglund, 
'Gay [94] 

SG 

Kronzon [94]~ 
KR. . 

Carter [13] 
- CA 

1 ftresen-t ' 
Study. 

in? 

+*/?• 

3-6 

.3 - 20 

.07 - ,17 

,3 - .6 

2 


,15 

.2 -7 

.07 

.05 - .1 

.125 


20 - 40 

1-2 

1-2 

3-7 

16 

Gridpoints 

1000 - 2000 

2000 -3500 

2000 -7000 

2000 - 5000 

2000 

Accuracy j 

A/ io" 6 

~ io" 6 

io"°- io" 6 

io” 5 

io"°- io" 4 

Max, numb! 
if. Iterats. 

5000 -9000 

2000 - 8000 

500 - 1000 

1500 -3000 

200 


layer 


UR, SG, MC - interaction of shock wave entering the boundary 

KR — interaction of flow due to a rear step 

CA - interaction of flow due to a forward compression corner 


NOTES : 


tl m *r fchod Mc relates to a grid that is fine and close to 
the boundary layer (while in the coarse, outer grid Ay/6o 
= 1-5 - 3) • 

b. The wide range of change in Ax and Ay, in method SG, stems 

rom a change in the grid so that its mesh size becomes a 
minimum close to the plate center, but increases in both x 
and y directions as one moves toward the periphery. 

c. The accuracy criterion for convergence is not well enough 
defined to apply to the same variables and in the same way 

or the various methods, so that the comparison must remain 
a qualitative one only. 

The comparisons do bring out that the method used here re- 
quires a relatively shorter time for convergence (by a whole 
order of magnitude) with lower accuracy than in previous studies 
(.out, as will be shown, with no significant computation errors). 


116 




n 


Appendix 0* - Formula for Computation of Suction Required 
to Prevent Separation 


A test was made to find, from the results of runs for 
suction under various flow conditions, a connection between the 
suction required to prevent separation and the data for the 
separated zone. 


For that purpose the following parameters are defined (all 
non-dimensional) : 


Area of separated zone (the bubble) 


. . / 

«• » 
\ 
i 




(for depth of unit) 


(0-1) 


Function of the flow at the bubble focus (the most negative value) 




when defined as 


i m / “ \ 

* f I . . • , % 


(0-2) 


Suction velocity (required to eliminate the separated zone) 


(when suction is normal to the wall 


V-^ s J K,r,Af 


t' 


u'i 


(0-3) 


n 


Now we will compute the expression C, which is defined by 


for various flow conditions 








(0-4 ) 


Note even when ^ have dimensions the value of C will 

be non-dimensional because the dimensions of its components are 

{w.7 - •*/. ; [v] * ;[? ] * *s/»* 


From the computation table (on the following page) it 
becomes clear that, generally, the value for c is in the region 
C = 0.30 - 0.42 (except in one instance where C = 0.24. In 
that case Re is exceptionally low and it becomes difficult to 
locate the boundaries of the separated region accuractely, 
because of the small differences in the velocity component 
values ) . 
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r\ 


n 


Since the cases examined cover a wide range of flow 
conditions, it may be assumed that the value expresses a 
singular relation between the suction required to prevent 
separation and the parameters of the separated region. 

To sum it up: the suction required to prevent separation 

can be designed by means of this formula: 

Vvj..' 5 C. [Set £ •• % ( o- 5 ) 

COMPUTATION TABLE: 



Original flow conditions without suction 

{Suction 1 


1 



Bubble dimens 

, Mesh dimens 

, Parameters 



Flow parameters 

(accdg.to nr. 

in the grid 

for the bubble 

Jer Fig. 4 

Run 




of meshes) 



focus 

t 

No. 

Ho 

• 1 o 

o 



£X-10 4 

ay-10 5 


(I'll-. : 

Uo 

102 

2.0 

1.4 

3xl0 5 

11.5 

3.5 

8.00 

5.00 

-.0011 

,6511 

-.023 

103 

2.0 

1.9 

3xl0 5 

12.5 

4,0 

8.00 

5.00 

-.0015 

.72211 

-.029 

W.->> 

104 ' 

•210 

3. i.5 

3x10* 

J {%6' 

4.5 

' i. oo 

' ‘ 5.00 

-.0026 

! ' V9472 

:: -7035 ' 

111 

3.0 

3.15 

3xI0 5 

12.0 

5. 5 

8.74 

5.46 

-.0011 

.4794 

-.027 

112 

-•rh' :vA- 

hll2 

4,0. 

• t. •£•:•»*" « 

4.0 

3.15 

"■.v’-r.-v.v 

1.9 ■ 

3xl0 5 

3x10° 

10. 5 
8.5' 

.3.5. 

VJV'Jl'. a 4 . 

2.5 

9. 54 

iVi, j-c.- 

"9.54 

r ^96. 

5.96 

-.0005 

k r -- 

-7 0003 

. .3430. 
“'.3516 

..-.018 . 
-.012 

113 

4.5 

3.15 

3x10 D 

10.0 

3.5 

£.92 

6.20 

-j. 0004 

.3123 

-.016 

122 . j 

- 2.Q 

.1.4 . 

5x1 0 4 

6.0 

2.5 

11.40 

7.10 

-.0011 

.. 7342 

. -.015 ... 

123 j 

2.0 

1.4 

lxlO B 

12.5 

3.5 


3.33 

-.0010 

1 

.64641 

-.02-3 

124 i 

2.0 

1.4 

5xl0 5 

13.0 

4,0 

2.75 

1.72 

-.0007] 

. 5421 1 

-.029 j 



y.^ y.,.V 

Vi? XI o' 
2 

tA*' 

U, 

v-tCfyl. 

r 

(x v-o.oiU 

WU, 

102 

6.44X10 -4 


-.023 

1.69 

3 * 45 



2.73 



3.96 

103 

8,00 


-.029 

2.08 


104 

If .30 

. 3.29 

-.035 

2.74 

4.20 

111 

12. SO 

. 3.55 

-.027 

2.29 

4.17 

112 

8 36 

2.C9 

-.013 

1.44 

3.65 

«112 

4. S3 

2.20 

-.012 

0.85 

3.09 

113 

e ei 

2.93 

-.016 

1.28 

3.63 ' 

122 

4,35 

2.21 

-.015 

1.50 

2.04 

123 

3,10 

1 .76 

-.026 

1.55 

r 






2.95 

124 
.. 

0.9S 

0.99 

-.029 

1.09 

» 

3.03 
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Appendix P - Description of the Computer Program 


a. Description of the Basic Program 

Pig. 4 shows a basic block diagram in which the various 
subroutines and their operation in the program are described. 

Below we will describe the main program and the subroutines 
attached to it. 

1) Main Program MAIN 

In this program control of the computer run procedure is 
exercised through the calling up of various subroutines. In it 
there are also definitions for a number of parameters that de- 
termine the desired options for computation, the form of the 
output and also the possibility for loading or unloading results 
into/from the assembler. Details of the parameters are given 
in paragraph c. 

2) ARVIS - Artificial Viscosity 

This subroutine is called up from GRID and its function is 
to add terms for artificial viscosity to the expressions in the 
difference schemes (in both integration steps). The coefficients 
of artificial viscosity, ccx and ccy, are determined by MAIN. 

3) CNBC - Boundary Condition at the Entry Section 

This subroutine computes the conditions at the entry section, 
which remain fixed in the course of the iterations and which are 
only dependent on the flow data Mo, Re , x , p and T . The 

X g S x O 

computation is carried out per the Polhausen method. 

4) DATA - Processing of Flow Data and Grid Preparation 

This subroutine calls up flow data like Mo, Re , x , p 

x s s 

and p /p and also the definitions of the computational field 
dimeniions through ratios !• .^Vr . . • From these 

data the dimensions of the computation grid are calculated, as 
well as the time interval for integration (based on the criterion 
for steadiness). 

5) DISC - Loading/Unloading of Intermediate Results into/from 
Assembler 

Through designation of Nip and/or NI for the I/O units this 
subroutine is employed to write and/or call up the appropriate 
results from the assembler on the disc. 



6) PLOW - Accompanying Calculations from the Results 
Obtained 

In this subroutine the Mach number and the flow function 
for the whole field, as well as the friction and heat transfer 
coefficients along the plate, are computed from the results of 
the computation. 

7) GRID - Operation of the Finite Difference Schemes 

In this subroutine the two-stage integration for the 
Brailovskaya difference scheme with second order finite dif- 
ferences is carried out, in which the variables u, v, p, t 
are calculated. (The subroutine includes a possibility for 
omitting the section along the shock wave in the case of the 
non-continuous computation method). 

8) INCND - Computation of Initial Conditions 

This subroutine computes the initial conditions for all 
the variables in the entire field and this is done based on 
considerations that were determined for the choice of initial 
conditions . 

9) OTPT - Printout of Results 

Since the printout of the results is a required condition 
this subroutine produces results for the variables in the compu- 
tational field x, u, v, p,p,t, y, m, ip 

(The friction and heat transfer coefficients are printed by the 
FLOW subroutine). 

10) PP - Pressure Computation 

This function computes the pressure from the equation of 
state. 

11) SHOCK - Computation of the External Boundary Condition 
Next to the Shock Wave 

This subroutine computes the boundary conditions on both 
sides of the shock wave, at the external boundary of the flow, 
and also permits calculation of the boundary conditions on 
both sides of the shock wave, along its entire length, by means 
of the non-continuous method in which the shock wave is removed 
from the computational field. 

12) STAB - Computation of Steadiness Coefficients 

In this subroutine the steadiness coefficients for each 
iteration are computed and they are printed out in tables, ac- 
cording to existing printout conditions (see note). 
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13) VRBC - Computation of Boundary Conditions (Except 
for the Shock Wave) 

In this subroutine those boundary conditions are computed 
that require updating for each iteration. They are located: 
along the plate, in the downstream flow and in the external 
boundary (except for downstream of the shock wave). 

14) YY - Computation of Viscosity 

This function computes the viscosity from Sutherland's 
formula. 

Note: Because STAB requires a very large memory and is not 

absolutely necessary for each computation, it can be dispensed 
with to save a lot of room and computer time (the method is 
described in paragraph c). 

b. Description of Additional Options 

1) The Difference Scheme by MacCormack 

The subroutine GRIDM (which exists parallel to GRID) 
performes the finite difference computation according to Mac- 
Cormack's method. (there is no possibility for the use of 
artificial viscosity and, therefore, no way to remove the shock 
wave from the computational field as is done in GRID). 

2) Difference Computation to the Fourth Order of Accuracy 

The subroutine GRID1 (which exists parallel to GRID) 
permits computation of the finite differences in the x-direction 
below the shock wave to the fourth order of accuracy, in 
contrast to the other parts of the field where second order 
accuracy prevails. 

3) Use of a Grid with Non-Uniform Meshes in the x-Direction 

The computation program can also be used for a grid with 
non-uniform meshes in the x-direction (which decrease in the 
direction towards the center) when Ax min =1.6Ay 

in the center of the grid. For that purpose the following 
subroutines exist (parallel to the subroutines of the uniform 
grid) : 

Subroutine DATAX instead of .DATA to process flow data and 
grid preparation. 

Subroutine GRIDX instead of GRID to compute finite differences 
(per Brailovskaya method only). 
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Subroutine VRBCX instead of VRBC to compute the boundary 
conditions . 

In addition to these subroutines there are two more functions 
HH and HL, by means of which the location and size of the 
meshes in the x-direction are determined. 

Note : The two functions HH and HL must be included in the 

program for the use of a uniform grid as well, since there are 
subroutines common to both types of grids (like CNBC, SHOCK 
and others) that require the use of those functions. 

c. Operation of the Program 


1) Control Card 

A batch type run of the program includes the following 
control cards: 

// JOB ... 

// EXEC . . . 

// FORT .SYSIN DD x 
MAIN program cards 
/x 

data cards 

// 

Parameters of the JOB card: 

Memory: without the use of STAB (250-350) K 

with the use of STAB (500-700) K 

Time for 200 iterations - without the use of STAB (200-300) sec 

with the use of STAB (400-500) sec 

Lines of print for 200 iterations including the printing of the 
initial conditions - 

without the use of STAB about 1700 lines 
with the use of STAB about 2300 lines 
(the ranges for change in the parameters are for the use of 
various options). 

For that reason it is recommended to use the following 
values with the JOB-card: 

without the use of STAB - T60, L02, R350 
with the use of STAB - T 85 , L03, R700 



2) Preparation of MAIN for a Run 


There are a number of indices in the MAIN program that 
control the procedure of the run and they must be defined 
before each run, in accordance with the requirements: 

NI - Use of a data assembler as initial condition: 

NI = 1 without use of assembler initial conditions are 
computed with INCND. 

NI / 1 initial conditions are read from assembler 
(whose logic number is identical to NI). 

NO - Use of data assembler to store output results: 

NO = 1 without use of assembler to conserve output 

NO ^ 1 results of the final iteration IMAX are loaded 
into the assembler (whose logic number is iden- 
tical to NO); the condition NO ^ NI must also 
exist . 

JAV - Use of artificial viscosity: 

JAV = 0 without use of artificial viscosity 
JAV = 1 with use of artificial viscosity 


OCX 

- Coefficients of artifical viscosity (generally in 
CCY the range 0.01 - 0.1) 


MODE - Use of continuous and non-continuous computation of 
the shock wave path: 

MODE = 1 continuous computation (regular computation 
method) 

MODE = 2 non-continuous computation (does not produce 
good results at this time). 

KPRT - Use of uniform and non-uniform grid: 

KPRT = 0 uniform grid 

KPRT = 1 non-uniform grid 

Q - Ratio of adjacent mesh lengths in a non-uniform grid 
(generally recommended Q = 1.1) 





IPRO - Possibility for printout of initial condition: 

IPRO = 0 without printout of initial condition 

IPRO = 1 with printout of initial condition 

I - Number of the first iteration in the run: 

1=0 when the run starts with initial conditions 
calculated through INCND 

I > 0 when the run starts with previous conditions 
called up from the assembler 

NC - Approximate boundary of the sonic line: 

This index indicates the value of the n-th point in the 
y-direction, in which the approximated sonic line exists (for 
computations of the non-continuous method) 

SSH - Initial value for the pressure ratio on both sides 
of the incident wave (this value is to a first approximation equal 
to the root of the overall pressure ratio between the downstream 
and upstream flow). 

Note : For the basic computation method it is required that 

KPRT = 0, MODE = 1, JAV = 1. Additional parameters that are 
not in use during this stage are introduced as fixed values 
into all runs. They are: 

SK = 1, J = 0, ISP = 0, INCD = 2, IGRD = 1, NCI = 2 

3) Preparation of Additional Data for DATA (or DATAX) 

In accordance with the input format the following data 
must be inputted: 

First card - 


EMO - Mach number of the external flow 

RE - Reynolds number at the intersection of the continuing 
shock wave with the plate 

PR - Prandtl number 

TPS - Initial value for the overall pressure ratio between 
downstream and upstream flow 

XSH - Distance, of the point of intersection between the 

continuing shock wave and the plate, from the leading 
edge 


/122 


124 



SY - Constant of the Sutherland equation (SY = 110°C) 




ALT - Reference altitude (M T ), used actually only to 

obtain the external temperature (In Its dimensional 
value in °C), which is required for the viscosity 
formula. (Generally ALT = 0 is used). 


Second card - 


G - The constant y = 1.4 

G1 - Ratio between the field length and the boundary layer 
thickness at entry section 

G2 - Ratio between the field width and the boundary layer 
thickness at entry section 

Z1 - Ratio between the grid mesh length and the boundary 
layer thickness at entry section 

Z2 - Ratio between the field width and the boundary layer 
thickness at entry section 

EC - Criterion for convergence (for the moment an exag- 
geterated value of EC = 1 x 10- 20 is introduced so that the 
program should stop before reaching IMAX) 

SD - Criterion for divergence (for the moment an exaggerated 
value of SD = 1 x 10 20 is introduced; usually, when there is a 
tendency for divergence the program will stop at a much earlier 
stage when negative values for temperatures are obtained, which 
causes calculation of the root of a negative expression in the 
viscosity formula, so conditions for termination of the run and 
printout of the results are inputted into GRID when the tempera- 
ture has a minus sign). 

Third card - 


IMAX - The maximum number of required iterations (generally 
no more than 200, so as to avoid overly long runs; when a larger 
number of iterations is required it is recommended to divide 
the run into stages, with the results of each stage loaded into 
the assembler and serving as initial conditions for the next 
stage ) . 

KPR - The number of iterations between printouts (recom- 
mended is KPR < 100). 


r ' 
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4) General Notes 

a) Before the run It should be verified that the ratios 
61/Z1 and 6 2/Z2 (whose Integral value is the number of grid- 
points in the x and y directions, respectively) match the size 
of the defined dimensions in COMMON and in DIMENSION. 

b) From experience it is known that the initial value 
of the pressure ratio between the downstream and the upstream 
flow, TPS, must be very little smaller (about 2 - 5 %) than the 
required final value. 

c) When the operation of subrouting STAB must be pre- 
vented (for reasons of saving space and computer time) a 
fictitious subroutine called STAB must be introduced together 
with MAIN, which will be designated: 

SUBROUTINE STAB 

COMMON/B/ 

IDIV = 0 

ICONV = 0 

RETURN 

END 

This subroutine (which will eliminate the original STAB) 
will generate the values ICONV, IDIV, which will prevent the 
needless stopping of the program. 

In place of this operation the callup CALL STAB can also 
be eliminated in MAIN and in its place to define IDIV = 0 and 
ICONV =0. As far as the computation is concerned there is 
nearly no difference between those two alternatives. 

d) To run the non-uniform grid (instead of the uniform 
grid) the following cards in 'MAIN must be changed: 

CALL DATAX instead of CALL DATA 

CALL GRID' instead of CALL GRID 

CALL VRBCX instead of CALL VRBC 

e) To use fourth order differences below the shock wave 
(instead of second order differences) the following card in 
MAIN must be changed: 



CALL GRID 1 instead of CALL GRID 

f) To use the MacCormack scheme instead of the 
Brailovskaya scheme, the following card must be changed in MAIN 

CALL GRIDM instead of CALL GRID 

Also, at the end of the subroutine DATA the command DT = DTY 
must be introduced. 
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TABLE 1A. ' CHECK OF COMPUTATION METHOD 


c 

Flow conditions for all runs: p /p = 1.4; Re x = 3*10 ; 


M =2.0 
o 


| Run No. 


ts-rA 


.j«y_ 



Notes 

rl • 

o 

o 

200 

0 

r — — — — 

l 

0 

0 

1 

Reference Run 

003 

200 

0 

2 

0 

0 

1 


005 

350 

0 

X 

i 

0 

1 

Artificial viscosity only for the 

■first. 1^0 Iterations, 

006 

200 

0 

1 

0 

1 

1 


007 

200 

0 

1 = 

0 

0 

2 

Component downstream of the shock \ 
wave's whole lenrth .computed hv l 

008 

200 

0 

2 

0 

0 

2 

average approximation saIne as a ^ ve 

009 

200 

0 

1 

0 

1 

2 

Isolation after 150 iterations 

021 

400 

1 

1 

0 

0 

1 


022 

400 

1 

1 

0 

1 

1 









Computation is carried out with 

031 

400 

1 

1 

0 

0 

1 

double accuracy (13 digits instead o:' 
f>-7 digits ) 

041 

.600 

1 

1 

0 

0 

1 

Computation ' carried out with fine grL 
14-7x4-8 (instead of the regular . 

051 

400 | 

1 

1 

0 

0 

1 

computation is carried out by Integra 




= - 

= = === 

= 




rs 


Notation: 


MODE 


= 1 continuous computation 
= 2 noncontinuous computation of the shock 
wave path 


KPRT = 0 unlform S r:i - d 

= 1 nonuniform grid in x direction 

TflV = 0 without artificial viscosity 
= 1 with artificial viscosity 


IGRD 


= 1 the entire field with second order differences 
= 2 part of the field below the shock wave with 
fourth order differences 


I S Tb = 0 without computations 'Of stability 
■ 1 with computations of stability 


IMAX - maximum number of iterations for the run 
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TABLE IB. CHECK OF INFLUENCE OF ARTIFICIAL VISCOSITY 


Run No. 

CC-- 

it— 

II 

II 

1 

II 

f. C / s 

:3sssssssa= 


Ho 

081 

.05 

.01 

3xl0 5 . 

1.4 

2 

082 

.10 

.01 

3x10 5 

1.4 

2 

083 

.10 

.05 

3xl0 5 

1.4 

2 


======== 

n 

ii 

ii 

n 

n 

ii 

n 


,=====, 

issas: 

091 

.10 

.01 

lxlO 6 

1.4 

2 

092 

.10 

.01 

5xl0 4 

1.4 

2 

. 093 

== bS 

1 

1 

O 1 

i 1 

• | 
U 

N 

.01 

3xl0 5 

3.15 

4 


Note: The first 150 iterations in this series of runs were 
performed with artificial viscosity, the following 200 without it. 




rs 
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TABLE 2. CHECK OF DEPENDENCE OF RESULTS ON PERIMETERS 

(All runs were performed according to the basic method 
•of run 001, up to 200 iterations) 


I Run Mo. 

101 

1 .72 

1 *s 'm 
.05 

3xlC S 

1.2 

Mo 

2.0 

Motes 

l> 

1 These results are used for 

102 ^ 

.72 

,05 

3xl0 5 

1.4 

2.0 

L comparison with previous 
/ analytical and experimental 

103 

.72 

.05 

3x10" 

1.9 

2.0 

\ results 

104 

.72 

.05 

3x10" 

_ e 

3.15 

2.0 


105 

.72 

.05 

3x10 

3.7 

2.0 

; No solution available for the 
"boundary conditions of the chara- 
cteristics since a subsonic. re/lio 
is dQ.\msireai!i_of the shock wave 

111 

.72 

.05 

3xl0 5 

3.15 

3.0 


112 

.72 

.05 

3x10 5 

3.15 

4.0 


113 

.72 

.05 

3x10 5 

3.15 

4.5 


114 

.72 

.05 

3x10 5 

' 

3.15 

5.0 

The solution is isolated after 
150 iterations 

121 

= s = = 

.72 

.05 

ixio 4 

1.4 

2.0 


122 

.72 

.05 

5x1 0 4 

1.4 

2.0 


123 

.72 

.05 

lxlO 6 

1.4 

2.0 


124 

.72 

.05 

5xl0 6 

1.4 

2.0 

From a practical point of view 
the flow is already turbulent and 
_the ^results therefore questionable 

131 . 

.72 

,03 

1.8X10 5 

1.4 

2.0 


132 

.72 

.07 

4.2X10 5 

1.4 

2.0 


133 

.72 

.09 

5.4xlO S 

1.4 

2.0 

1 

141 

l.C 

= s ss =3 e = 

.05 

i ii 

i n 

1 Ol II 

1 X II 

1 H* II 

l O II 

i cn ii 

i ii 

i ii 

> ■ H 

1.4 

2.0 

:ssse: 



*Run 102 is similar to 001. 

o 


/l4l 
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TABLE 3- INFLUENCE OF BOUNDARY CONDITIONS (THERMAL CONDUCTIVITY 
AND MASS TRANSFER BY THE PLATE) ON THE FIELD OF FLOW 
’ (All runs were performed according to the basic method 
of run 001, up to 200 iterations). Flow conditions 
for all runs: P e /p Q = 1.4; Re Xg = 3-105; M q = 2.0. 




O 



Run No. 

acj aaaitsa 
<201 

1 S uc tion region w. ref. to start 



1 End 3 eirt 

Lnnlri'T 

V-t 

J)ist, 

- a. ruu 

6,0 

Grid pt. 

: s*j.1 

is 

GO 

Gist. 

3.2 

Grid pt. 

: r a * ss a a a 
\ aauJiaa a 

15 

: a = a ai a 
0 

-0.03 

UJ* 

issaaia 

0.4 

202 

0,6 

203 

0.0 

204 

1.0 

205 

1.2 

106 

211 

1.4 

<cJT r ) 

;o 

,v'ss' 

1 

1 

212 

6,0 

GO 

3.. 2 

15 

-0.02 

213 

C.3 

. GO . 

3.2 

15 

-0.01 

. 214 

G.O 

CO 

3.2 

15 

0.01 

215 

221 

C.G 

r.iajc. 

6.0 

-60 

;‘J3SJauJ 

GO 

3.2 

r aa ac a r 
2.6 

15 

ataiKUwa a 

0 

0.02 

-0.03 

222 

C.O 

GO 

2.0 

0 

-0.02 


223 

G.O 

GO 

2. G 

0 

-0.01 


224 

G.O 

GO 

3.0 

22 

-0.03 


OOr. 

J 

6.0 

GO 

3.0 

22 

-0.02 


226 

G.O 

GO 

3.0 

22 

-0.C1 


227 

C.3 

CO 

4.4 

30 

-0.03 

i 

220 

6.0 

CO 

4.4 

30 

-0.02 

■ 

220 

G.O 

a r j «* _ 

GO 

<i a J .1 U J S i 

4.4 

r u u l- &> 

30 

a a u w m » u i: 

-0.01 
n b a :i • u l 


N otes 


j5-iiaaaa4 J3« 


Check of influence of 
heating or cooling of 
plate surface on the 
interaction 


iJ l 


Check of influence of 
suction or injection at 
the plate surface on 
the interaction 


vxsii&juissaKxaceiaaa 


Check of influence of 
location for start of 
suction on the size 
of the separated regioni 


i 

i 

I 

i 

i 
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TABLE 4. INFLUENCE OF LENGTH AND LOCATION OF SUCTION ZONE ON THE /l43 
SEPARATION 

•(All runs were performed according to the basic compu- 
tation method of run 001, up to 200 iterations). Flow 

conditions for all runs: Re Y = 3«105; n /n = 14 . 

A s ’ *0 5 

M 0 = 2.0. (The length of the separated region without 

suction, under these conditions, is in the range of 
35 < m < 45 grid points on the plate surface.) 



f Suction 

Region w. ref. to the Start f 

Hun i<o® 

! ' r- nU . . 

0 uarr g 


Distance 

. ; ' !-..i . 

Grid Point 

Distance 

Grid Point 

251 

4.0 

25 

3.2 

15 

252 

4.8 

35 

3.2 

15 

253 

5.6 

45 

3.2 

15 

254 

6.4 

55 

3.2 

15 

255 

4.8 

35 

4.0 

25 

256 

5.2 

40 

4.0 

25 

257 

5.6 

45 

4.0 

25 

258 

6.0 

50 

4.0 

25 

259 

4.8 

35 

4.4 

30 

260 

5.2 

40 

4,4 

30 

261 

5.6 

45 

4.4 

30 

262 

6.0 

50 

4.4 

30 

263 

5.2 

40 

4.8 

35 

264 

5.6 

45 

4.8 

35 

265 

6.0 

50 

4.8 

35 

266 

5.6 

45 

5.2 

40 

267 

: : = = SS s =. S =; 

6.0 
= = = = = 

50 

5.2 

ll 

II 

II 

II 

O II 

^ It 

II 
II 
II 
II 
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TABLE 5A. INFLUENCE OF SUCTION ON THE SEPARATION UNDER VARIOUS /144 
FLOW CONDITIONS 

(All runs were performed according to the basic compu- 
tation method of run 001, for 200 iterations.) A 
region of suction is present in the range of 15 < m < 

< 60 grid points on the plate surface for all runs. 


Rim No. 

VU/ 

do 

X s 


Me 


t | 

; Notes 1 

301 

-0.03 


3xl0 5 

2.0 

1.9 


302 

-0.02 


3xl0 5 

2.0 

1.9 


303 

-0.01 

.05 

3xl0 5 

2.0 

1.9 

■I m Li6n c 6 ox duc wioii on x>xic 

304 

-0.03 


3xl0 5 

2.0 

3.15 

-Lnuensiby oi vaxious ohock ft aves 

305 

-0.C2 


3X10 5 

2.0 

3.15 

1 

306 

-0.01 . 


3xl0 5 

2.0 

3.15 


311 

-0.03 


5x1 0 4 

2.0 

1.4 


A 312 | 

: 

-0.02 


5x1 0 4 

2.0 

1.4 

Influence of Suction on the 
various Reynolds Numbers 

313 

-0.01 

.05 

5x1 0 4 

2.0 

1.4 

314 

-0.03 

lxlO 6 

2.0 

1.4 

315 

-0.02 


lxlO 6 

2.0 

1.4 


n 

ii 

ii 

ii 

II w 
II M 
» O) 

II 

l 

-0.01 


lxlO 6 

2.0 

1.4 

ii 

ii 

ii 

ii 

n 

n 

ii 

ii 

ii 

i 

ii 
ii 
u 
ii 

i 

ii 
ii 

i 

i 

i 

i 

i 

i 

i 

i 

i . 
i 
i 
i 

i , 
i 
i 
i 

i 

321 

-0.03 


3xl0 5 

3.0 

3.15 


322 

-0,02 


3xl0 5 

3.0 

3.15 

Influence of Suction on the 

323 

-0.01 

.05 

3xl0 5 

3.0 

3.15 

324 • 

-0.03 


3x10° 

4.0 

3.15 

various Mach Numbers 

325 j 

-0.02 


3xl0 5 

4.0 

3.15 


326 | 

-0.01 

l 

1 

1 

I 

1 

1 

3xl0 5 

4.0 

3.15 

1 

11 

II 

II 

11 

II 

II 

II 

II 

It 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

II 

If 

II 

1! 

II 

II 

II 

II 

11 

II 

II 


n 


148 
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TABLE 5A (continued) PARAMETRIC CHECK FOR VARIOUS 


Run 

No. 


/ X * 
C 


Ho 

P */Po 1 

Notes *. 










331 

-0.005 


3xl0 5 

2.0 

1.4 

i 

Influence of suction on the 









332 

-0.015 


3x10° 

2.0 

1.4 

various shock wave intensities 


333 

-0.005 


3xl0 5 

2.0 

1.9 





0.03 

c: 





334 

-0.015 


3x10 

2.0 

1.9 



335 

-0. 005 


3xl0 5 

2.0 

3.15 

» 

'\ 


336 

-0.015 


3xl0 5 

2.0 

3.15 



341 

-0.005 


5x1 0 4 

2.0 

1.4 

Influence of suction on the 


342 

-0.015 


,5xl0 4 

2.0 

1.4 

various Reynolds numbers 




0.03 






343 

-0.005 


1x10° 

2.0 

1.4 



344 

-0.015 

o 

lxlO 6 

2.0 

1.4 

' * 


351 

-0.005 


3xl0 5 

3.0 

3.15 

\ 


352 

-0.015 


3xl0 5 

3.0 

3.15 

\lnfluence of suction on the 

■K 

0.03 

cr 



various Kach numbers 


353 

-0.005 


3x1 (j 

4.0 

3.15 



354 

-0.015 


3xl0 5 

4.0 

3.15 

• 

sssrts 

s csaSQ 



3x10^ 





361 

-0.015 


2.0 

1.4 



362 

-0.025 





Influence of suction on the 



3x10 

2.0 

1.4 

various shock wave intensities 


363 

-0.015 


3xl0 5 

2.0 

1.9 



364 

-0.025 

0.07 

3x10 5 

2.0 

1.9 



365 

-0.015 


3xl0 5 

2.0 

3.15 



366 

-0.025 


3xl0 5 

2.0 

3.15 

• 


371 

-0.015 


5x1 0 4 

2.0 

1.4 

1 

I 


372 



4 



Influence of suction on the j 
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TABLE 5B. CHECK OF VARIOUS COMBINATIONS OF FLOW CONDITIONS 
DURING SUCTION 

• (All runs were performed according to the basic compu 
tation method of run 001, up to 200 iterations). The 
suction zone is located for all runs in the range 
15 < m < 60 grid points on the plate surface. 
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Fig. 2 - Boundary Conditions 
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Fig. 11 - Distribution of Stability 
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Fi£. 15 - Conrparlson of Run 1Q1 to nVW 
Experimental and Comput- 
ational results 
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Fig* 23 - Hach nunber H. influenco 
oa tho Interaction. 
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Fig. 2h - Flowfield for Run 1 1 2 112 
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Fig. 28 - Prandtl number ? r influence 
on the Interaction. 
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Fig. 33 - Influence of variation 
of oTw/u. 1 • on the 

Boundary Layer shape. 
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Fig. 36 - Dependence of separation length ^*on suction 
velocity for various Mach numbers 
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